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COMPUTER  MODELS  FOR  TWO-DIMENSIONAL 
TRANSIENT  HEAT  CONDUCTION 


Mary  Remley  Albert 


INTRODUCTION 

Most  major  Army  installations  are  heated  with  central  heat  distribution  systems.  The  thermal 
regimes  around  buried  distribution  systems  are  of  interest  for  several  reasons,  not  the  least  of  which 
is  the  estimated  millions  of  dollars  lost  annually  from  these  systems  (Phetteplace  et  al.  1981). 

Also,  when  replacing  damaged  insulation  or  installing  insulation  in  a  new  system,  it  is  desirable  to 
know  what  an  optimum  balance  is  between  the  initial  cost  of  insulation  and  the  continued  oper¬ 
ating  cost  of  heat  losses.  In  addition,  freezing  and  thawing  of  t^e  ground  around  a  buried  distri¬ 
bution  system  could  damage  it,  through  loss  of  support  and  settlement,  for  example. 

Analysis  of  the  thermal  regime  involves  the  solution  of  the  heat  conduction  equation  for  situ¬ 
ations  with  complicated  geometries  and  a  variety  of  boundary  conditions.  Usually,  these  partial 
differential  equations  cannot  be  solved  analytically;  we  must  resort  to  numerical  methods.  One 
long-established  method  is  that  of  finite  differences,  a  relatively  straightforward  numerical  method 
used  successfully  to  solve  a  variety  of  differential  equations.  A  finite  difference  computer  program, 
set  up  in  a  general  form  to  solve  the  heat  conduction  equation  under  a  variety  of  geometries  and 
boundary  conditions,  provides  a  powerful  tool  with  which  the  engineer  can  assess  problems  in 
conductive  heat  transfer. 

The  objective  of  this  paper  is  to  document  the  development  and  verification  of  two  general 
two-dimensional  finite  difference  computer  programs  that  were  written  to  model  time-varying 
heat  conduction  in  a  medium  composed  of  many  materials.  The  first  program.  ADI,  solves  for 
the  general  case  in  which  no  change  of  state  occurs  in  the  conducting  medium.  The  second  pro¬ 
gram,  ADBPC,  is  an  adaptation  of  ADI  that  includes  the  effects  of  phase  change  in  the  conducting 
medium.  The  programs  are  able  to  handle  convective,  constant  flux,  specified  temperature  and 
semi-infinite  boundaries.  Material  properties  such  as  thermal  conductivity,  density  and  specific 
heat  are  allowed  to  vary  with  time  or  temperature.  The  programs  are  relatively  easy  to  use,  that 
is,  they  are  easily  set  up  for  new  conduction  problems  by  those  who  have  a  background  knowledge 
of  FORTRAN. 

The  computer  programs  were  written  and  programmed  in  FORTRAN  by  the  author  on  CRREL’s 
PRIME  400  computer. 


FINITE  DIFFERENCES  APPLIED  TO  HEAT  TRANSFER 


Finite  differences  are  commonly  used  in  the  numerical  solution  of  partial  differential  equations. 
The  procedure  involves  the  replacement  of  differentials  by  differences,  and  is  best  illustrated  by 
construction  of  the  so-called  finite  difference  grid.  Let  us  examine  a  grid  that  is  set  up  in  two  di¬ 
mensions  to  model  heat  conduction  (Fig.  1).  The  two  dimensions  in  this  case  represent  spatial 
independent  variables,  x  and  y.  Each  point  of  the  grid  is  called  a  node  and  represents  the  area 
enclosed  by  the  square  around  it  (with  unit  depth).  T(x,y)  represents  the  temperature  of  node 
x,y.  For  each  node,  the  temperature  and  material  properties  are  assumed  uniform  for  the  region 
it  represents.  By  specifying  the  initial  temperatures  at  each  of  the  nodes,  the  nodal  thermal  proper¬ 
ties  and  the  boundary  conditions,  we  can  solve  the  heat  conduction  equation  to  determine  the 
temperatures  of  the  nodes  at  later  times. 


r 


- 9 - 

- — • - 

-  9 - 

> 

• 

• 

TU,y-Ay) 

• 

• 

• 

• 

• 

• 

T(x-Axty! 

• 

T(x.y) 

• 

T(x*Ax,y) 

• 

• 

• 

• 

• 

• 

T(x,y+Ay) 

• 

• 

• 

• 

• 

• 

• 

• 

• 

* 

• 

< 

Ay 


Ax 


Figure  1.  A  simple  finite  difference  grid. 


Heat  conduction  equation 

The  equation  governing  transient  heat  conduction  in  two  dimensions  is 


dT 

dr 


(1) 


where 


x  and  r  =  spatial  variables 
T-  temperature 
t  =  time 

k  =  thermal  conductivity  (this  may  vary  over  the  spatial  and  time  domains) 
p  =  density 
C  =  specific  heat. 


The  nonlincarities  in  this  partial  differential  equation  that  arise  from  the  inclusion  of  the  phase 
change  condition  will  be  discussed  in  the  Phase  Change  section. 

Partial  differential  equations  may  be  expressed  in  finite  difference  form  cither  by  a  Taylor’s 
series  expansion  about  a  point  or  by  physical  considerations,  such  as  a  heat  flow  balance  in  the 
case  of  the  heat  conduction  problem.  For  problems  involving  variable  thermal  properties  or  com¬ 
plicated  boundary  conditions,  the  heat  balance  approach  is  the  simplest. 
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Figure  2.  A  node  on  the  interior  of  the  grid. 


Consider  the  control  volume  represented  by  a  node  on  the  interior  of  the  grid,  illustrated  in 
Figure  2.  The  heat  flow  Q  into  the  control  volume  from  any  adjacent  node  may  be  calculated  as 
follows: 

0=  l  A  AT  (2) 

where  k  -  thermal  conductivity 

L  =  distance  over  which  the  heat  flow  occurs 
A  =  area  (per  unit  depth)  perpendicular  to  the  direction  of  heat  flow 
A T  -  change  in  temperature  between  the  two  nodes. 

Examine  the  heat  flow  from  node  N  to  node  P.  The  conductance  k/L  is  the  reciprocal  of  the  re¬ 
sistance  between  the  two  nodes.  The  resistance  between  the  two  nodes  RNf  is  simply  the  sum  of 
the  individual  resistances, 

^NP  =  ^N  +  *p  ' 

Each  resistance  is  the  distance  divided  by  the  material  conductivity, 

R  =M2  +  M2 

NP  *N  *p 


*  =  -J_  =  -2-*p*-p—  (3) 

L  R^p  +  kp) 

In  the  computer  program  it  is  assumed  that,  for  nodes  not  on  a  boundary  of  the  grid,  A^  and  Ax 
are  uniform  and  equal  throughout  the  grid.  Let  the  nodal  spacing  Ax  =  Ay  =  As.  Following  the 
form  of  eq  2,  we  find  that  the  heat  flows  into  the  region  associated  with  node  P  from  nodes  N, 
W,  S  and  E  may  be  calculated  as  follows. 


(4) 


2k  k 

^  =  A^/tJ+fcp)  (As)(7w  -  ^p)  (5) 


Csp  = 


2kskp 
As(fcs  +  kp) 


(Ai)(rs-rp) 


(6) 


Qe p  = 


Ikpkp 
Aj(Arg  +  kp ) 


(As)(7’E-7’p) 


(7) 


Note  that  the  area  perpendicular  to  the  direction  of  heat  flow  is  taken  per  unit  depth,  that  is,  A  - 
As  •  1. 

The  sum  of  the  heat  flows  into  a  node  is  responsible  for  the  temperature  change  of  the  node: 

Qnp  +  £?wp  +  (?sp  +  @EP  =  ^p  ‘  aJ  (ftCpTp-pCpTp)  (8) 

where  Vp  is  the  control  volume  for  unit  depth  (  Vp  =  (As)2  in  this  case)  and  the  primed  variables 
represent  time  t  +  At. 

Now  substitute  eq  4-7  into  eq  8  to  get  the  finite  difference  form  of  eq  1  for  a  node  not  on  the 
boundary  of  the  grid: 

^NP^N  "  T’p)  +  kyjp  (Tw  -  Tp)  +  &Sp  (T%  -  Tp)  +  kE(7'E  -  Tp)  =  ^  (p  CpTp  -  p  CpTp). 

(9) 


Since  As  is  uniform  for  interior  nodes  of  the  grid,  substitutions  like  the  following  have  been  made 
for  simplicity’s  sake: 


k 


NP  " 


*N  +*P 


(10) 


Equation  8  describes  the  temperature  change  for  one  node  P  over  one  time  step.  The  values  p, 
Cp  and  k  are  held  constant  for  each  time  step,  but  may  be  changed  between  steps.  It  is  necessary 
to  apply  eq  8  (or  the  appropriate  boundary  equation)  to  each  node  of  the  grid  to  solve  for  the 
temperature  distribution  for  one  time  step.  The  process  is  repeated  until  the  desired  number  of 
time  steps  have  been  completed.  Let  us  now  investigate  several  methods  used  in  the  solution  of 
the  set  of  equations. 

The  first  and  probably  the  easiest  method  to  program  is  the  explicit  method,  where  eq  8  is 
solved  for  Tp  (which  is  the  value  of  Tp  for  time  t  +  At).  Note  that  all  of  the  other  temperature 
variables  represent  time  t,  so  that  the  equation  may  be  solved  explicitly.  An  initial  temperature  is 
assigned  to  each  node  in  the  grid;  the  appropriate  equation  is  solved  for  each  node  to  determine  the 
temperature  resulting  from  the  first  time  step.  The  resultant  temperature  distribution  is  then  used 
in  calculations  for  the  second  time  step.  This  process  continues  until  the  desired  number  of  time 
steps  have  been  calculated.  The  main  problem  with  this  method  is  that  a  stability  criterion  must 
be  met  to  produce  a  reliable  answer.  For  a  two-dimensional  homogeneous  grid,  the  criterion  is 


Note  that  a  balance  must  be  struck  between  the  internodal  distance  (As)  and  the  time  step  (At). 

In  problems  involving  steep  temperature  gradients,  a  small  internodal  spacing  must  be  used;  con¬ 
sequently,  the  time  step  must  be  very  small.  The  combination  of  many  nodes  and  time  steps  makes 
the  solution  costly  in  computer  time. 

A  second  method  used  in  solving  finite  difference  equations  is  the  implicit  method.  Here  Tp 
is  taken  as  the  value  of  Tp  for  time  t,  and  all  of  the  other  variables  represent  time  t  +  A f;  hence, 
all  temperature  variables  are  unknown  except  Tp.  An  appropriate  equation  must  be  written  for 
each  node  in  the  grid,  and  the  resulting  set  of  equations  is  solved  simultaneously  for  each  time  step. 
For  the  two-dimensional  case  of  an  x  by  y  grid,  there  will  be  xy  equations  in  xy  unknowns.  The 
matrix  of  coefficients  is  banded  (bandwidth  =  2x  +  1 ),  but  in  general  there  is  no  symmetry  about 
the  main  diagonal.  Large  matrices  often  occur  and  they  are  usually  solved  by  an  iterative  method. 
The  implicit  method  has  the  advantage  of  being  unconditionally  stable,  but  it  may  require  a  fair 
number  of  iterations  for  adequate  convergence  in  the  matrix  solution. 

The  finite  difference  method  chosen  for  this  computer  program  is  known  as  the  alternating 
direction  implicit  method  (Peaceman  and  Rachford  1955,  Carnahan  et  al.  1969).  It  is  uncondi¬ 
tionally  stable  like  the  fully  implicit  method,  yet  it  requires  only  the  solution  of  a  tridiagonal 
matrix,  thus  being  efficient  in  computer  time  and  computer  core  storage.  (The  matrix-solving  al¬ 
gorithm  is  discussed  in  the  TRIDIG,  The  Matrix  Solver  section).  Essentially,  the  method  requires 
the  solution  of  two  different  equations  for  each  time  step.  The  first  equation  is  implicit  only  in 
the  horizontal  direction.  The  results  from  the  first  step  are  then  used  to  solve  the  second  equation, 
which  is  implicit  only  in  the  vertical  direction. 

For  example,  apply  the  alternating  direction  considerations  to  a  node  in  the  interior  of  the  grid. 
Equation  9,  which  describes  the  heat  balance  for  an  interior  node,  is  repeated  here  for  convenience: 


*np  (T’n  "  ^p)  +  *wp  (rw  "  ^p)  +  ^sp  (7-s  “  ^p)  +  ^ep  (^e  "  ^p) 

(A  sYp'Cp  ^  (As)2  pCp 

At  Tr  '  At  Tv 

Each  node  in  the  grid  is  given  an  initial  temperature.  Then,  for  the  first  pass,  let  the  above  equa¬ 
tion  be  implicit  in  the  horizontal  direction.  The  time  step  for  the  pass  will  be  At/2,  and  the  above 
equation  may  be  written  as  follows: 


*wp?w 


+  *EP^E 


|fcwp  +fcEP  + 


2(As)2  p'Cp 


~*NP^N  "  *SP^S  + 


2(As)2pCp  1 

*np+*sp”  £t  \  • 


(11) 


The  primed  variables  represent  the  value  of  those  variables  at  time  t  At/2.  The  right-hand  side 
of  the  equation  is  known;  7^,,  Tp  and  T'e  are  to  be  determined.  Equation  1 1  (or  the  appropriate 
boundary  equation)  is  applied  to  each  node  of  a  row  of  the  grid,  creating  a  tridiagonal  matrix  of 
coefficients.  Such  a  system  of  equations  is  quickly  solved  for  each  row  of  the  grid.  The  resulting 
temperature  distribution  represents  t  +  At/2,  the  end  of  the  first  pass. 

The  process  is  similar  for  the  second  pass,  except  that  now  eq  9  is  rewritten  to  be  implicit  in 
the  vertical  direction : 


*np^n  +  *sp^s  ~ 


*np  +  *sp  + 


2(As)2p'Cp 


5 


The  primed  variables  represent  time  t  +  At,  and  the  others  represent  time  t  +  At/2.  The  tridiagon¬ 
al  matrix  is  formed  and  solved  for  each  column  of  the  grid;  the  resultant  temperatures  represent 
the  temperature  distribution  in  the  grid  after  t  +  At,  one  time  step.  The  entire  process  is  repeated 
until  the  desired  number  of  steps  have  been  calculated. 

The  reader  may  consult  Carnahan  et  al.  (1969),  Holman  (1972),  Croft  and  Lilley  (1977)  and 
Mitchell  and  Griffiths  (1980)  for  more  information  on  finite  difference  methods. 

Boundary  conditions 

The  following  boundary  conditions  will  be  derived  from  heat  balance  considerations.  In  this 
paper  each  boundary  equation  will  not  be  expressed  in  the  form  needed  for  use  in  the  alternating 
direction  procedure,  but  the  reader  should  be  aware  that  each  equation  was  put  into  that  form  for 
use  in  the  computer  programs  ADI  and  ADIPC. 

Sides  of  the  grid 

Consider  a  node,  P,  on  the  right-hand  grid  boundary.  The  control  volume  associated  with  the 
node  is  the  area  enclosed  by  dotted  lines  around  it,  as  depicted  in  Figure  3. 

The  sum  of  the  heat  flowing  through  the  boundaries  of  P's  control  volume  is  responsible  for 
the  temperature  change  of  the  node,  that  is. 

Ax  Ay  pCp  AT? 

Qnp  +  Q\mp  +  Qsp  +  Q  ip=  (13) 

where  <2NP  =  Ihe  heat  flow  from  node  N  to  (or  from)  node  P 
Ax 

2  ■  Ay  =  the  volu  ne  (for  unit  depth)  of  node  P 
p  =  the  density  of  the  material  in  node  P 
C p  =  the  specific  heat 
T  =  temperature 
t  =  time. 


The  flow  of  heat  from  node  N  to  node  P  is  given  by 


As 


Figure  3.  A  node  on  a  right-hand  boundary. 


where  *NP  is  the  effective  conductivity  between  nodes  N  and  P,  and  for  Ax  =  Ay  =  As, 

,  2*n*p 

NP~  (*n+*p) 

Similarly,  the  heat  flow  from  nodes  W  and  S  to  node  P  may  be  given  as 

Gwp  =  *wp(7w  '  rp)  (*5) 

2s p  =  J  *sp  (Ts  -  Tp)  .  ( J  6) 

Now  consider  several  cases  describing  the  heat  flow  across  the  boundary. 

Constant  flux  boundary.  For  a  boundary  subject  to  a  constant  heat  flux,  the  condition 

bT  I 

aT  =  constant 
dx  lx=P 

exists  at  the  boundary.  For  node  P,  illustrated  in  Figure  3,  the  heat  flow  across  its  east  side  is 
given  by 

QEr  =  <t>Ay  (17) 

where  <f>  is  the  heat  flux  per  unit  area  crossing  the  boundary.  The  equation  for  a  node  on  the  right- 
hand  side  of  the  grid  with  constant  flux  is  obtained  by  combining  eq  13-17  and  allowing  Ax  = 

A_y  =  As: 

J  *np  (T’n  -  TP)  +  *WP  (7W  -  TP)  +  kSP  (Ts  -  rp)  +  *As  =  ^  [p'  c;t;-P  CpTv ] . 

(18) 

The  indices  may  be  suitably  rearranged  for  constant  flux  boundaries  on  other  sides  of  the  grid. 

Note  that  for  a  boundary  which  is  insulated  or  on  a  line  of  symmetry  the  zero  heat  flux  condi¬ 
tion  holds,  and  0  =  0. 

Convection  boundary.  For  the  node  illustrated  in  Figure  3,  exposed  to  convection  on  the  right- 
hand  side,  the  heat  flow  from  the  convective  medium  to  node  P  is  given  by 

Qv.v  ~  hp.&.v  (Fe  -  Tp)  (19) 

where  hE  is  the  coefficient  of  convective  heat  transfer  and  TE  is  the  temperature  outside  the  grid. 
Combining  equations  13,  14,  15,  16  and  19,  we  arrive  at  the  equation  for  a  node  on  the  convective 
right-hand  side  boundary,  for  Ax  =  Ay  =  As: 

T  ^np  (^n  ”  ^p) +  ^'wp  (T’w  ”  ^p)  +  ~2  *sp  (T's  '  Tp)  +  hAs  (TE-Tp)  = 

^  Ip’CpT;  -pCpTp  ].  (20) 

Specified  temperature  boundary.  For  a  node  of  specified  temperature  on  a  boundary,  corner 
or  inside  tire  grid,  apply  the  equation  Tp  =  C,  where  C  is  the  temperature  of  the  node  at  time  t. 

Semi-infinite  boundary.  This  condition  represents  a  continuous,  uniform  material  extending 
in  one  direction,  with  a  known  temperature  a  large  distance  away.  It  is  approximated  here  by  use 
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Figure  4.  A  node  on  a  semi-infinite  boundary. 

of  a  large  intemodal  distance  between  the  last  two  nodes  of  a  row  of  the  grid.  Consider  the  situa¬ 
tion  illustrated  in  Figure  4,  for  the  right-hand  side  of  the  grid.  Node  E  is  not  actually  a  part  of 
the  grid,  but  is  the  location  of  known  temperature  TE  outside  the  grid. 

The  distance  between  the  last  two  nodes  and  also  between  the  boundary  node  and  the  location 
of  the  known  temperature  TE  is  ZljAs.  Note  that  represents  a  multiple  of  As.  The  heat  flow 
from  nodes  N,  W,  S  and  E  to  node  P  may  be  calculated  as  follows 

Gjyp  =  kjyp  (2Di  -  1)  ( Tn  -  Tp)  (21) 

£?wp  =  ^wp  (Tw  ~  Tp)  (22) 

C?sp  =  *sp  "  1)  (^s  -  rp)  @3) 

Cep  =  ^ep  (Te  ~  TP)  .  (24) 

The  effective  conductivities  are  of  the  same  form  as  eq  10,  except  for  kwp  and  kEP 

*WP=  [*p  +  (2Z>i-l)*w]  • 

Node  E  is  assumed  to  be  the  same  material  as  node  P,  thus  kEP  =  A:p.  The  sum  of  the  heat  flows 
accounts  for  the  temperature  change  of  the  node,  thus  the  equation  for  a  node  on  the  right-hand 
semi-infinite  boundary  is 

^NP  (2Di  ~  0 (Tn  -  Tp)  +  kyjp  (rw  -  Tp)  +  kSp  (2Dj  -  1)  (7g  -  Tv)  + 


*EP  (^E  "  Tp)  = 


(.2D,  -  1 )  (As)2 


(fi'CpTp  -pCp  Tp) 


where  (2£>j  -  1)  (As)2  is  the  control  volume  for  the  node  per  unit  depth. 

In  any  finite  difference  formulation,  the  accuracy  of  the  solution  increases  as  the  area  repre¬ 
sented  by  a  node  decreases.  Therefore,  when  using  the  semi-infinite  boundary  formulation,  the 
user  should  specify  the  smallest  D{  acceptable.  The  semi-infinite  condition  should  only  be  used 
in  regions  where  the  temperature  gradient  is  small  and  precise  knowledge  of  the  temperature  dis¬ 
tribution  is  not  critical. 
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Figure  5.  An  interior  node  adjacent  to  the  semi-infinite 
boundary. 


This  boundary  conditions  also  requires  a  special  heat  balance  for  the  node  adjacent  to  the  semi¬ 
infinite  node.  Consider  again  the  situation  for  a  right-hand  boundary,  as  illustrated  in  Figure  5. 
The  appropriate  heat  flow  equations  are 


Cnp  =  *np  V*  "  7p) 

QWp  =  kyfp  (Ty,  -  Tr) 

(?sp  =  ^sp  Vs  “  ^p) 

where  fcNP,  kwp,  ksp  are  of  the  same  form  aseq  10,  and 
Cep  =  *ep  Vn  ~  TP) 


'Lkpkp 

*«■  [*,♦(»,- d*,i  • 

The  equation  for  an  interior  node  adjacent  to  a  right-hand  semi-infinite  boundary  node  is 
^NP  ^p)  +  kyfp  (Tw  -  Tp)  +  k§p  (Tg  -  Tp)  +  fcEP  (Tg  -  Tp)  = 

^  [P‘c;t; -pCpTp] .  (3o> 

Corner  of  the  grid 

The  corners  of  the  grid  require  special  heat  balances,  depending  on  the  conditions  on  each  edge 
of  the  grid.  Consider  the  upper  right-hand  corner  of  the  grid,  shown  in  Figure  6.  T*  and  T?  are 
not  a  part  of  the  grid,  but  are  temperatures  outside  the  grid.  The  heat  flow  from  the  two  nodes 
adjacent  to  node  P  may  be  given  as 

Gwp  =  J  *wp  (^w  "  7p)  (31) 


CSp  =  2  *sp  Vs  ~  Tf) 


where  ky,p  and  kSP  follow  the  form  given  in  eq  10.  £NP  and  QPP  are  dependent  upon  the 
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Figure  6.  A  node  on  a  corner  of  the  grid. 


particular  boundary  conditions  and  will  be  presented  shortly.  The  comer  equations  will  follow  the 
form 

(As)2  AT 

C?np  +  @wp  +  Gsp  +  @ep  =  4  P  ^p  y  (33) 

(As)2 

where  is  the  volume  represented  by  node  P  for  unit  depth. 

Constant  flux  on  both  sides.  Let  ^  be  the  heat  flux  per  unit  area  crossing  the  north  side  of  the 
comer  shown  in  Figure  6,  and  0E  be  that  crossing  the  east  side  of  the  corner.  Then  the  heat  flows 
across  the  two  sides  are  given  by 

C?n  p  =  0n  '  y  (34) 

£?ep  =  0e  ’  y  ■  (35) 

The  equation  for  the  corner  is 

2  *wp  (7-w  "  7p)  +  J  *sp  <rs  “  ^p)  +  J  +  =  ^4Ar  ^  ^p  ^p  ~  p  ^p  ^pl 


As  in  the  case  of  the  side  of  the  grid,  if  a  side  of  the  corner  is  insulated  or  on  a  line  of  symmetry, 

0  =  0  for  that  side. 

Convection  on  both  sides.  For  a  corner  subject  to  convection  on  both  of  its  sides,  the  heat  flow 
from  each  of  the  two  sides  may  be  given  as 

Onp  =  an  y  (^n  ~  ^p)  (^7) 

0EP  =  *e  y  (Te-Tp)  (38) 

where  hN  is  the  coefficient  of  convective  heat  transfer  on  the  north  side  and  hE  is  that  on  the  east 
side.  The  equation  for  node  P  is 

2  *wp  (Tw  ~  Tf)  +  -j  kSP  ( Ts  -  Tp)  +  —  hNAs  (TN  -  Tp)  +  —  hEAs  (TE  -  Tp)  = 

^  [ p'C’T ;  -PCpTp)  .  (39) 


»>>  .> .*• 
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Figure  8.  A  semi-infinite  side  node  adjacent  to  a  semi¬ 
infinite  comer. 


Except.for  this  change  in  fcNP,  the  equation  for  the  square  interior  node  inside  an  upper  right-hand 
comer  with  semi-infinite  conditions  on  both  sides  is  the  same  as  eq  29. 

Now  consider  the  semi-infinite  nodes  adjacent  to  the  semi-infinite  corner  node.  These  are  node 
temperatures  Ty,  and  Ts  in  Figure  7.  The  figure  is  presented  again,  this  time  as  Figure  8,  with  the 
nodes  relabeled.  The  type  of  node  under  consideration  is  represented  by  nodes  A  and  P  in  Figure 
8.  The  following  equations  reflect  the  heat  flow  to  node  P: 


0NP  *NP  ~  0  (^N  ~  ^p) 

(45) 

CwP  =  *WP  (7-w  "  ^p) 

(46) 

08P-*8P  (20,-1)  (7,-7,) 

(47) 

Cep  =  *ep  (Ts  ~  ^p)  • 

(48) 

Here, 

2kNkp 

1.  —  n  r  -  -  J  #.  — 

2kyfkp 

*NP  *,,+(2^-1)*,,  and*WP=  kp  +  ilD;-  l)*w  ■ 


The  equation  for  this  node  is 

p  -  1 )  (  -  7"P)  +  kw  p  (7"w  -  TpJ  +  fcgp  (2Z?j  -  1 )  (7"g  -  7"p)  + 

(2D,  -  1)  (As)2 

*E p  (Tr  -  Tp)  =  - - |p'  C'  rp  -  P  Cp  Tp]  .  (49) 

Constant  flux,  convective  corner.  For  the  corner  illustrated  in  Figure  6,  allow  convection  to 
occur  across  the  north  side  of  node  P  and  constant  flux  across  the  east  side 


Gnp  =  *n  y  (Tti  ~  Tf)  (50) 

QE p  is  given  by  eq  35,  and  Qwp  and  QSP  are  given  by  eq  31  and  32,  respectively.  The  resulting 
equation  for  the  corner  is: 

2  (T'n  "  ^p)  +  J  Vw  ~  T?)  +  2  ^sp  (Ts  ~  ^p)  +  J  = 

[p'c;rP-PcfTf] .  (so 

Constant  flux,  semi-infinite  comer.  The  corner  illustrated  in  Figure  9  has  a  semi-infinite  bound 
ary  on  the  right  side.  Let  there  be  a  constant  flux  per  unit  area,  0N,  across  the  north  side  of  node 
P,  then 


(?NP  ~  '  0 

Gwp  =  J  *wp  ^ 

Gsp  =  ^sp  (^Dj  ”  0  (Ts  ~  Tp) 
Qep  =  2  ^EP  -  ^p^ 
*"p“  *,  +  (20,-1)** 


it 


sp  = 


2kskP 
k «,  +  kp 


^ep  ~  *p 


(52) 

(53) 

(54) 


(55) 


The  equation  for  the  corner  is  as  follows: 

(2Dj  -  1 )  As  +  (Tw-Tp)  +  kSP(2Di-l)(Ts-TP)+~  (TE-TP)- 


(2D-\)(As)2 

— jAt -  [P'c;t;.pcptp]. 


(56) 


Figure  9.  A  node  on  a  semi-infinite  corner. 
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Again,  the  node  adjacent  to  the  semi-infinite  node  must  have  a  special  equation.  For  node  W  in 
Figure  9,  with  constant  flux  ^  across  the  top,  eq  57  applies,  with  referance  to  the  node  under 
consideration  as  node  P: 

As  +  j  ky, p  (Tw  -  Tp )  +  kSP  ( Ts  -  Tf)  +  j  kEp  (TE  -  Tp)  = 

\\f  Ip  CpTp  -  p  Cp  Tp\  .  (57) 

i 

s 

Semi-infinite,  convective  comer.  For  the  comer  in  Figure  9,  allow  convection  across  the  north 
side,  and  allow  the  same  semi-infinite  boundary  on  the  east  side, 

Gnps*n  (20|  -  1)(TN  “  Tp)  •  (58) 

Owp.  GSP  and  <2ep  appear  in  eq  53, 54  and  55,  respectively.  Thus  the  equation  for  the  configura¬ 
tion  is 

kv,  p  k  pp 

hN  (ZD,  -  1)  (rN  -  Tp)  +  -Y  (Ty,  -  Tp)  +  *SP  (2D,  -  1)  (Ts  -  Tp)  f  -f-  (TE  -  Tr)  = 


(2D,  - 1)  (As)2 
2  At 


Ip'c;t ;  -pCpTp] . 


(59) 


The  appropriate  equation  for  node  W  in  Figure  9,  subject  to  convection  across  the  top,  follows. 
Again,  refer  to  this  node  as  node  P: 


An  As  (Tn  -  Tp)  + 


2  ^wp  (T’w  ~  ^p)  +  *SP  Vs  ~  ^p)  + 


Vy,  ~  TP)  = 


2  Ar 


[p’C^T;  -p  Cp  Tp]  . 


(60) 


PHASE  CHANGE 

The  methods  discussed  in  this  report  so  far  apply  to  both  AD1PC  and  ADI,  the  programs  de¬ 
veloped  for  heat  conduction  with  and  without  phase  change.  The  ideas  on  phase  change  presented 
in  this  section  will  apply  only  to  ADIPC.  The  computer  program  solves  the  heat  conduction  equa¬ 
tion  only;  possible  effects  of  moisture  migration  through  the  medium  (and  unfrozen  water  content 
in  a  frozen  soil,  for  example)  are  neglected.  Although  the  examples  used  involve  freezing,  ADIPC 
permits  either  freezing  or  thawing  at  any  node  of  the  grid. 

During  the  freezing  process,  the  temperature  of  a  substance  that  is  initially  above  its  fusion 
temperature  decreases  as  heat  is  removed  until  the  substance  reaches  its  fusion  temperature.  The 
continuing  removal  of  heat  produces  no  change  in  temperature  until  the  energy  equivalent  to  the 
latent  heat  of  the  substance  is  removed.  The  substance  has  then  changed  state  to  become  frozen; 
any  further  removal  of  heat  again  results  in  a  decrease  of  temperature  below  the  temperature  of 
fusion. 

ADIPC  models  freezing  or  thawing  by  defining  an  apparent  specific  heat  during  phase  change 
that  accounts  for  the  entire  enthalpy  change  that  takes  place,  including  the  enthalpy  in  the  latent 
heat  of  fusion.  To  do  this,  it  must  be  assumed  that  phase  change  takes  place  over  a  finite  tempera¬ 
ture  range,  AT,  around  the  fusion  temperature.  Allow  Tf  to  denote  the  fusion  temperature,  Cp 


and  **8  to  be  the  specific  heat  of  the  substance  above  and  below  freezing,  respectively,  and  let 
Cf*  be  the  apparent  specific  heat  defined  for  phase  change.  For  a  temperature  change  in  a  subsi 
not  involving  a  change  of  state,  the  enthalpy  change  A H  is  given  as 


A// 


-/ 


CPdT. 


Now  the  temperature  range  for  phase  change  will  be  Tt  ±  A  772,  and  the  apparent  specific  heat  ac¬ 
counts  for  the  entire  enthalpy  change  as  follows,  where  Ht  represents  the  latent  heat  of  fusion: 


i 


T  +  Ar 

Tf  T 


AT 


Ct  dT = 


7>-t 


1 


r  AT  rB 
7f"  2 


T  +  AT 
ft  * 

Cp_  dT*  I  Cp  dT  +  7/l 


/: 


(61) 


A7'is  chosen  small  enough  so  that  each  specific  heat  may  be  assumed  constant  over  the  integral, 

A  T  AT 

C?47-=CPb  -  *CrA  t  *Hl 

Cp*-i(CPBtCPA)t^.  (62) 

Equation  62  defines  the  apparent  heat  capacity  as  used  in  ADIPC.  The  apparent  heat  capacity 
method  is  further  discussed  in  Bonacina  and  Comini  (1973). 

Phase  change  is  implemented  in  ADIPC  as  follows.  At  the  beginning  of  each  time  step,  each 
node  of  the  grid  is  examined.  If  its  temperature  lies  within  the  range  Tf  ±  A  772,  the  specific  heat 
for  that  node  is  defined  by  eq  62.  The  conduction  equations  are  set  up  and  solved  as  usual.  At 
the  end  of  the  complete  time  step,  the  newly  calculated  temperature  of  each  node  is  compared 
with  the  temperature  of  that  node  at  the  beginning  of  the  time  step.  If  the  temperature  skipped 
from  below  Tf  -  AT  17  to  above  Tf  +  A 772,  or  vise  versa,  then  the  phase  change  front  skipped  that 
node,  and  the  temperature  of  the  node  is  reassigned  as  follows.  If  the  temperature  of  the  node 
skipped  from  the  frozen  to  the  unfrozen  domain, 


where  T'  is  the  calculated  nodal  temperature  before  its  reassignment.  If  the  node  skipped  from  an 
unfrozen  state  to  the  frozen  state,  its  temperature  is  reassigned  as  follows: 


In  this  way,  the  program  assures  that  the  phase  change  front  does  not  skip  a  node. 

A  drawback  with  the  apparent  heat  capacity  method  is  that  it  is  not  designed  to  follow  the 
exact  location  of  the  phase  change  front  for  each  time  step;  it  is  designed  to  calculate  nodal  tem¬ 
perature.  The  location  of  the  T(  isotherm  may  be  found  by  interpolation,  but  it  approximates 
the  location  of  the  front  in  a  step-like  pattern  because  of  the  discretization  of  a  continuous  space, 
as  will  be  illustrated  in  the  program  verifications.  This  behavior  is  minimized  by  use  of  a  small 
internodal  spacing.  The  method  is  flexible  enough  to  deal  with  phase  change  in  two-dimensional 
space  when  there  may  be  several  locations  of  phase  change  fronts;  for  this  flexibility  the  method 
was  chosen. 


COMPUTER  PROGRAM 


In  using  finite  differences,  we  replace  differentials  and  derivatives  with  differences.  It  makes 
sense,  then,  that  the  accuracy  of  the  solution  increases  as  these  differences  become  smaller.  Spe¬ 
cifically,  the  accuracy  improves  as  the  nodal  spacing  and  time  step,  As  and  At,  are  made  small. 
There  is  no  set  criterion  for  just  how  small  they  must  be,  usually  this  is  discovered  through  trial 
and  error. 

The  grid  for  these  programs  may  be  any  size,  so  far  as  the  outside  dimensions  are  concerned, 
but,  except  for  the  semi-infinite  boundary  conditions,  each  node  in  the  grid  represents  a  square 
area  of  unit  depth.  RAY  (I,  J,  K)  is  the  array  that  represents  the  grid.  I  and  J  are  spatial  variables 
in  Cartesian  coordinates;  I  is  incremented  vertically  and  J  horizontally.  Khas  three  values.  RAY 
(I,  J,  1 )  represents  the  temperature  distribution  in  the  grid,  i.e.,  a  temperature  is  stored  at  each  I, 

J  location.  RAY  (I,  J,  2)  records  the  nodal  location  type.  Examples  of  location  types  include  a 
node  on  a  constant  flux  boundary,  variable  interior  node,  specified  temperature  corner  node,  etc. 
Each  location  type  is  assigned  a  number,  and  one  such  number  is  stored  for  each  I,  J  location  under 
RAY  (I,  J.  2).  This  information  is  used  to  assign  the  correct  form  of  the  heat  conduction  equation 
to  each  node.  Each  material  type  in  the  problem  is  given  a  number;  these  numbers  are  stored  in 
RAY  (I,  J,  3)  and  may  be  used  to  assign  a  conductivity,  density  and  specific  heat  for  each  node. 

The  programs  are  made  up  of  four  parts:  1 )  a  data-gathering  subroutine,  2)  the  main  program, 
3)  a  subroutine  to  solve  the  tridiagonal  matrix  and  4)  a  subroutine  to  locate  the  user-specified 
isotherms  at  times  specified  by  the  user.  Parts  3  and  4  are  identical  for  both  ADI  and  ADIPC,  and 
part  1  is  similar  for  each.  Therefore,  parts  1 , 3  and  4  will  be  discussed  only  once  but  are  included 
in  each  of  the  two  programs. 

ADDATA,  the  data  subroutine 

The  necessary  data  for  the  program  are  handled  mainly  through  subroutine  ADDATA.  Vari¬ 
ables  used  are  defined  in  the  comment  statements  at  the  beginning  of  the  subroutine.  The  user 
has  simply  to  edit  the  subroutine,  following  directions  in  the  comment  statements  in  the  subrou¬ 
tine  to  initialize  variables  and  arrays.  When  the  subroutine  is  run,  the  data  are  put  into  a  formatted 
data  file  ADIDAT;  the  user  does  not  have  to  worry  about  formats. 

TRIDIG,  the  matrix  solver 

This  subroutine  solves  the  tridiagonal  matrix  formed  in  the  main  program.  The  matrix  is  formed 
by  the  application  of  eq  1 1  or  12  (or  a  boundary  condition  counterpart)  to  each  node  of  the  grid. 
The  resultant  system  of  equations  is  illustrated  in  matrix  form  in  Figure  10.  The  matrix  of  coef¬ 
ficients  is  an  n  by  n  matrix;  all  entries  are  zero  except  those  on  the  three  center  diagonals,  hence 
the  term  tridiagonal  matrix.  The  vector  on  the  right  contains  the  elements  of  the  right-hand  side 
of  each  equation.  To  conserve  computer  storage  space,  only  the  three  diagonals  are  stored,  as 
vectors.  Thus  the  effective  storage  space  required  is  reduced  from  n  by  n  +  I  to  n  by  4. 

The  solution  algorithm  is  commonly  used  (e.g.  Gerald  1980).  Each  element  of  the  lower  diag¬ 
onal  may  be  eliminated  by  subtracting  the  appropriate  multiple  of  the  (i  -  1)  row  from  the  r'th 
row.  The  values  of  b}  and  dx,  after  elimination  of  <Zj,  are 

(65) 
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Figure  10.  Didiagonal  matrix  equation. 


After  replacing  the  element  of  vectors  b  and  d  by  the  new  values,  we  perform  a  back  substitution 
as  follows: 

<■-£  (67) 


"i’-T 


,  for  i  =  n  -  1 ,  n  -2, ...  1  . 


The  elements  of  the  solution  vector  replace  vector  d. 

ISOTHM,  die  isotherm  finder 

Subroutine  ISOTHM  examines  the  temperature  distribution  in  the  grid  and  performs  a  linear 
interpolation  between  adjacent  nodes  to  produce  the  Cartesian  coordinates  of  the  locations  of  the 
user-specified  isotherms.  The  coordinates  are  listed  in  file  POINTS  in  two  columns,  one  for  the 
horizontal  coordinate  and  one  for  the  vertical  coordinate.  POINTS  may  be  used  as  a  data  file  for 
a  plotting  routine.  Subroutine  ISOTHM  may  be  called  at  the  completion  of  any  number  of  time 
steps;  the  frequency  is  specified  by  the  user  in  ADDATA. 

ADI,  main  program 

ADI  first  interactively  asks  the  user  if  ADDATA  should  be  called.  The  subroutine  should  be 
called  if  ADIDAT  does  not  already  contain  the  formatted  data.  If  the  user  so  desires  ADDATA 
is  called.  If  not,  ADI  reads  the  data  from  ADIDAT. 

For  each  time  step,  the  thermal  properties  such  as  conductivity,  density  and  specific  heat  are 
updated,  if  appropriate,  and  the  resultant  conductivities  between  the  nodes  are  figured.  Each  node 
is  assigned  the  coefficients  for  the  appropriate  form  of  the  conduction  equation ,  forming  the  tri- 
diagonal  matrix  (Fig.  10),  and  subroutine  TRIDIG  is  called  to  solve  the  matrix.  The  resultant 
temperature  distribution  is  then  used  in  the  second  half-time  step,  and  TRIDIG  is  called  once 
again.  These  temperatures  then  represent  the  distribution  for  one  complete  time  step.  If  appro¬ 
priate,  subroutine  ISOTHM  is  then  called  to  locate  the  user-specified  isotherms.  This  procedure 
is  repeated  for  the  desired  number  of  time  steps. 

The  initial  data  and  boundary  conditions  are  printed  into  file  ADIOUT ;  the  temperature  dis¬ 
tributions  for  specified  time  steps  are  printed  into  file  ADITMP,  and  finally,  a  new  data  file, 
ADNDAT,  is  created  for  the  final  temperature  distribution.  This  file  is  in  the  same  format  as  the 
input  file  ADIDAT;  ADNDAT  may  be  used  as  a  starting  point  if  the  user  wishes  to  run  the  model 
for  more  time. 


^ 


ADIPC,  main  program 

ADIPC  is  similar  to  ADI,  but  differs  in  that  phase  change  considerations  are  implemented.  At 
the  beginning  of  each  time  step,  if  the  temperature  of  a  node  lies  in  the  phase  change  temperature 
range,  the  program  assigns  the  apparent  specific  heat  to  that  node,  as  previously  discussed.  Other¬ 
wise,  the  specific  heat,  thermal  conductivity  and  density  for  each  node  are  updated  according  to 
the  user’s  specifications. 

At  the  end  of  the  time  step,  the  temperature  of  each  node  is  compared  to  its  temperature  at  the 
end  of  the  previous  time  step.  If  the  temperature  skipped  over  the  phase  change  temperature  range, 
it  is  reassigned  as  discussed  in  the  Phase  Change  section,  "  hen,  if  appropriate,  subroutine  ISOTHM 
is  called. 

The  initial  and  boundary  conditions  are  printed  in  readable  form  in  file  ADPOUT;  the  tempera¬ 
ture  distributions  for  specified  time  steps  are  printed  into  file  ADPTMP,  and  a  new  data  file  with 
the  final  temperature  distributions  is  printed  into  file  ADPNDT. 

VERIFICATION  OF  ADI 
Comparison  of  ADI  with  analytical  results 

Semi-infinite  corner 

The  results  of  ADI  will  first  be  compared  to  the  problem  of  a  semi-infinite  corner,  as  illustrated 
in  Figure  1 1 . 


Figure  1 1.  Semi-infinite  corner  problem. 


Tt  is  the  uniform  initial  temperature;  at  time  /  =  0  the  temperature  of  the  two  edges,  given  by 
x  =  0  and  r  =  0.  is  changed  to  T0.  The  solution  is  well  documented  (Carslaw  and  Jaeger  1959,  and 
Holman  1972)  and  is  found  by  using  a  product  solution  for  two  one-dimensional  problems.  The 
solution  is  given  by 


(69) 


where  a  is  the  thermal  diffusivitv.  and 
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Figure  12.  Finite  difference  grid  for  the  semi-infinite  corner 
problem. 

For  a  comparison  run,  the  following  values  were  used: 

T.  =  20°  C  (68°  F) 

7'0=40°C(J04°F) 
a  =  0.0025  m2/hr 

The  finite  difference  grid  appears  as  shown  in  Figure  1 2.  The  boundary  condition  assigned  to  the 
top  and  left  side  was  that  of  constant  temperature  equal  to  40°C.  The  right  side  and  bottom  of 
the  grid  were  assigned  a  semi-infinite  boundary  condition,  with  Ts  =  Tf  =  20°C  and  Di  =  50  (see 
the  Sides  of  the  Grid  section).  The  nodes  not  situated  on  a  boundary  were  assigned  an  initial  tem¬ 
perature  of  20°C.  A  time  step  of  0.25  hr  was  used,  and  the  internodal  spacing  was  0.025  m. 

The  locations  of  the  35°,  30°  and  25°C  isotherms  are  plotted  for  the  results  of  ADI  and  for  the 
analytical  solution  for  several  time  steps  in  Figure  13.  Excellent  agreement  is  found  between  the 
two  solutions  for  regions  not  adjacent  to  the  semi-infinite  boundary.  As  previously  stated,  the 
semi-infinite  condition  is  an  approximation.  Also,  the  large  internodal  distance  used  for  the  semi¬ 
infinite  boundary  decreases  the  accuracy  of  the  solution  in  that  region. 


Figure  13.  locations  of  35°,  30°  and  25°  C  isotherms  from  ADI  and  from  the  analytical  solution  for  the  semi- 
infinite  corner  problem. 
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Figure  14.  Finite  rectangle  problem. 


Finite  rectangle 

In  order  to  provide  a  two-dimensional  comparison  that  doesn’t  involve  the  use  of  semi-infinite 
boundaries,  and  also  to  demonstrate  internodal  spacing  effects,  consider  the  problem  of  a  rectangle 
of  uniform  initial  temperature  T{  that  is  subject  to  a  step  change  in  the  temperatures  of  all  of  the 
edges  to  T0  at  time  zero.  The  problem  is  illustrated  in  Figure  14.  The  solution  of  this  problem  is 
again  obtainable  from  the  product  solution  of  two  one -dimensional  problems.  Carslaw  and  Jaeger 
(1959)  provide  the  one  dimensional  solution: 


T\x,t)-T0 

T,-T0 


i-ai 


The  product  solution  is  then  given  by 


T(x,y)~  T0  + 
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A  short  FORTRAN  program,  INFSUM,  was  written  to  compute  for  this  equation  the  various 
space  and  time  increments;  it  is  listed  with  the  output  in  Appendix  A.  For  the  comparison,  the 
rectangle  is  indicated  by  2  =  1  m,  b  =  2  m,  T0  =  40°C,  Tx  =  20°C,  and  the  thermal  diffusivity  a  is 
set  equal  to  0.00251  m2/hr. 

ADI  was  run  several  times  to  demonstrate  the  effect  of  internodal  spacing  on  the  accuracy  of  the 
solution.  Because  of  the  symmetry  of  the  problem,  we  have  to  model  only  one  quarter  of  the 
problem,  assigning  a  zero  flux  boundary  condition  to  edges  of  the  grid  that  fall  on  lines  of  sym¬ 
metry.  The  region  modeled  is  that  portion  of  the  rectangle  which  lies  in  the  fourth  quadrant  of 
the  Cartesian  graph.  The  first  run  was  made  with  an  internodal  spacing  of  0.1  m  and  a  time  step 
of  5  hr.  The  resulting  24°,  28°,  32°,  36°  and  40°C  isotherms  are  plotted  for  several  time  steps 
(Fig.  15a).  An  inspection  of  the  data  printed  in  ADITMP  reveals  that  the  maximum  discrepancy 
between  the  two  solutions  is  1 ,5°C,  and  it  occurs  in  the  first  time  step  (5  hr)  in  the  lower  right- 
hand  corner.  This  location  at  early  times  represents  the  steepest  temperature  gradient  in  the  prob¬ 
lem.  By  the  second  time  step  (10  hr),  the  discrepancy  reduces  to  a  maximum  of  0.6°C.  The  max¬ 
imum  difference  between  the  solutions  continues  to  decrease  as  the  temperature  gradient  decreases 
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b.  Internodal  spacing  0.05  m;  time  step  2.5  hr. 

Figure  15.  Locations  of  24°,  28°,  32°,  36°  and  40°C  isotherms  (temperature  increasing  from  left  to  right) 
from  ADI  (finite  difference)  and  analytical  solution. 

and  the  location  of  the  maximum  temperature  difference  changes  to  the  locations  of  the  steepest 
gradients.  The  input  and  output  for  this  run  of  ADI  are  listed  in  Appendix  A. 

The  same  problem  was  run  with  an  internodal  spacing  of  0.05  m  and  a  time  step  of  2.5  hours. 
The  resulting  isotherm  plots  are  shown  in  Figure  15b.  A  closer  agreement  was  found  between  the 
results  of  ADI  and  the  analytical  solution.  A  node-by-node  comparison  shows  the  maximum  dif¬ 
ference  to  be  1  °C  for  the  5-hr  distribution,  decreasing  to  0.4°  in  the  10  hr  distribution.  In  gen¬ 
eral,  the  accuracy  of  the  solution  is  improved  as  the  internodal  spacing  decreases;  the  user  must 
determine  the  accuracy  demanded. 


One-dimensional  semi-infinite  problem 

It  is  also  of  interest  to  compare  the  results  of  ADI  to  the  one-dimensional  problem  of  a  medium 
initially  at  a  uniform  temperature;  the  surface  temperature  then  undergoes  a  step  change  to  a 
different  temperature,  and  the  resulting  temperature  distribution  is  examined  over  time.  This 
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Figure  16.  One-dimensional  semi-in  finite  problem. 


problem  corresponds  to  the  theoretical  setup  of  the  experimental  problem  examined  in  the  Com¬ 
parison  of  ADI  with  Experimental  Results  section.  The  problem  is  illustrated  in  Figure  1 6. 

The  solution  to  this  analytical  problem  is  found  in  many  texts  on  heat  transfer;  the  reader  may 
refer  to  Holman  (1972).  The  solution  is  given  as 


T(x.t)-T0 

T,-T0 


=  erf 


x 

2  y/ar 


(72) 


where  x  =  the  distance  below  the  surface 
t  =  time 

a  =  the  thermal  diffusivity,  0.00251  m2/hr  in  this  example. 

The  grid  used  for  this  was  3  by  40  nodes,  with  an  internodal  spacing  of  0.025  m  and  a  time 
step  of  0.25  hr. 

For  two  depths,  0.05  and  0.20  ni,  the  results  of  this  analytical  solution  are  compared  to  the  re¬ 
sults  of  ADI  over  9  hr.  The  results  are  given  in  Table  1 . 


Table  1 .  Comparison  of  ADI 
results  to  analytical  results  for 
the  one-dimensional  semi-in¬ 
finite  problem. 


The  results  compare  to  within  0.02°C  for  times  after  2  hr.  Earlier  times  have  temperatures  that 
differ  by  as  much  as  0.24°C  near  the  top  of  the  grid.  Probably,  this  is  because  the  steepest  temper¬ 
ature  gradients  occur  in  this  problem  soon  after  the  step  change  in  temperature  occurs  at  the  sur¬ 
face.  Recall  that  steeper  temperature  gradients  require  a  smaller  intemodal  distance;  if  we  want 
greater  accuracy,  the  model  could  be  again  run  using  an  intemodal  spacing  of  less  than  0.025  m. 

The  smaller  the  intemodal  distance,  the  more  accurate  the  solution. 

Comparison  of  ADI  with  experimental  results 

As  an  example  of  the  use  of  ADI,  let  us  compare  it  to  some  experimental  results.  Data  are 
available  from  a  full-scale  experiment  that  was  conducted  in  a  4-  by  5-m  container  of  soil  that  was 
1 .2  m  deep.  Of  the  4-m  width,  only  the  center  1 .76  m  was  included  in  the  experiment;  the  rest 
was  outside  an  insulated  boundary  and  was  put  there  for  support.  The  soil  used  was  Lebanon  sand 
of  19.7%  moisture  content.  The  experiment  was  designed  to  start  at  an  initial  uniform  tempera¬ 
ture  of  70°F  (21.1 1°C)  in  the  sand  in  the  box,  then  the  surface  temperature  was  to  be  raised 
100°F  (37.78°C),  and  the  change  in  the  temperature  distribution  over  time  was  to  be  monitored. 

An  analytical  solution  is  available  for  this  case,  and  ADI  was  also  run  for  this  case.  The  results 
compare  very  well  (see  the  One-Dimensional  Semi-infinite  Problem  section,  analytical  comparisons). 

In  the  comparison  of  the  experimental  data  to  calculated  results,  temperatures  were  used  that 
represent  a  20-cm  wide  band  taken  vertically  through  the  center  of  the  box.  Three  strings  of 
thermocouples  were  placed  in  this  band,  as  shown  in  Figure  17. 


O  I  m 


Figure  1 7.  Thermocouple  locations  in  the  experimental  setup  ( top  of  soil 
container  on  right). 


Figure  18.  Surface  temperatures  from  experiment  used  in  ADI  runs. 


In  the  actual  experiment,  the  initial  temperatures  ranged  from  20.6°  to  22.4°C.  The  initial  tem¬ 
perature  distribution  for  ADI  was  taken  as  an  average  of  the  three  thermocouple  string  temperatures 
for  each  depth.  It  required  approximately  27  hr  to  raise  the  surface  temperature  to  37.78 °C,  after 
which  the  surface  temperature  remained  in  the  range  from  37.2°  to  38.4°C.  Because  the  data  logger 
malfunctioned,  12  hr  of  data  are  missing  in  a  period  beginning  9  hr  after  the  start  of  the  experiment; 
however,  the  rest  of  the  equipment  continued  to  operate  normally.  In  comparisons  made  with  these 
data  it  is  assumed  that  the  surface  temperature  increased  linearly  over  time  for  this  period.  See 
Figure  18  for  a  graph  of  the  surface  temperature  vs  time  for  the  data  and  that  used  in  the  compari¬ 
son  run  with  ADI. 

Soil  tests  were  conducted  at  CRREL  to  determine  the  density  and  thermal  conductivity  of  the 
sand.  The  density  was  determined  to  be  1996.25  kg/m3;  this  density  was  assumed  constant  in  the 
ADI  comparison.  The  conductivity  was  measured  at  two  temperatures  above  freezing.  At  4.44°C 
the  conductivity  was  determined  to  be  1 .673  W/m  K,  and  at  26.67°C  it  was  1 .803  W/m  K.  Each  node 
in  the  ADI  model  was  assigned  a  conductivity  according  to  its  temperature  (in  degrees  C)  for  each 
time  step  from  the  following  equation: 

kWi  =  0.005857'+  1.647  .  (73) 

This  equation  was  determined  from  the  two  conductivity  tests. 

The  value  of  the  specific  heat  was  taken  from  measurements  done  on  Lowell  sand  by  Kersten 
(1949).  Again,  each  node  in  the  ADI  model  was  assigned  a  temperature-dependent  specific  heat. 

The  equation  fit  to  Kerstcn's  data  is 

CP  =  0.0039r+  0.34927  .  (74) 

i.j 

The  2  by  48  grid  used  in  ADI  had  a  distance  between  nodes  of  0,025  m.  The  temperature  at  the 
surface  was  specified,  and  was  taken  from  the  data  previously  discussed.  The  time  step  was  0.25  hr; 
the  surface  temperature  for  each  time  step  was  interpolated  linearly  from  the  data.  The  sides  and 
bottom  of  the  grid  were  assigned  the  zero  heat  flux  boundary  condition.  The  problem  was  modeled 
for  48  hr. 

The  graphs  in  Figure  19  compare  the  results  for  several  times  (9  hr,  24  hr  and  48  hr)  during  the 
run;  they  show  a  reasonable  agreement  but  the  model  tended  to  underpredict  the  temperature 
change. 
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Figure  19.  Results  of  two  ADI  runs  using  different  conductivities  plotted  against  experimental  results. 

The  model  was  then  run  again,  using  conductivity  measurements  on  Lowell  sand  by  Kersten 
(1949).  The  temperature-dependent  equation  for  conductivity  used  for  this  run  was 

k. j-  0.005857-  +  1.975  .  (75) 

These  results  may  again  be  found  in  Figure  19.  This  time  excellent  agreement  was  found  between 
the  actual  data  and  the  ADI  calculations. 

Each  of  the  two  runs  of  ADI  required  2  minutes  and  1 3  seconds  of  computer  time  on  the  PRIME 
computer  at  CRREL. 

VERIFICATION  OF  ADIPC 

Comparison  of  ADIPC  with  analytical  results-the  Neumann  solution 

ADIPC  was  first  compared  to  the  well-known  one-dimensional  analytical  Neumann  solution.  A 
semi-infinite  region  is  initially  at  a  uniform  temperature  T0  which  is  above  the  fusion  temperature 
Tr  Suddenly,  the  surface  temperature  is  changed  to  Ts,  a  temperature  below  the  fusion  tempera¬ 
ture.  The  subsequent  movement  of  the  phase  change  front  may  be  calculated. 

Allow  p,  Cp,  k  and  HL  to  represent  the  density,  specific  heat,  thermal  conductivity  and  latent  heat 
of  fusion,  respectively.  Thermal  diff usivity  a  is  defined  as  k/pCp.  Temperature  7" is  a  function  of 
position  and  time.  Subscripts  I  and  2  refer  to  the  property  or  variable  in  the  frozen  and  liquid 
phase,  respectively.  The  boundary  conditions  follow: 


T ,  (0,r)  =  r5 

(76) 

Tj  ( x ,  o)  =  r0 

(77) 
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The  temperatures  in  the  liquid  and  solid  regions  must  satisfy  the  following  equations 

_l  9£,=0 

dx2  3* 


a2r2  j  ar2 

ax2  ’  *2 


The  solution  has  been  given  many  times  in  the  literature  (e.g.  Carslaw  and  Jaeger  1959).  The  solu¬ 
tion  for  the  location  of  the  phase  change  front  follows: 

X=2 X  s/aj  (82) 

where  X  must  be  determined  from  the  relation 
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For  comparison  the  following  values  were  used: 


**  =2'21  mK 
p,  =917  ^ 


=  0.5815 


*2  =0.580  ^ 


p2  =998.2 


CB  =1.16 


H ,  =93.00  -r— - 
L  kg 

7>  =  0°C  . 

Ts  =  -4.67°C 

T0  =  4.67°C 

The  problem  was  modeled  for  3  hr  using  an  internodal  spacing  of  0.5  cm  and  a  time  step  of  0.0025 
hr.  The  depth  of  the  0°C  isotherm  is  plotted  with  the  analytical  solution  in  Figure  20.  As  men¬ 
tioned  earlier,  the  location  of  the  phase  change  front  is  found  by  interpolating  between  the  nodal 
temperatures  to  find  the  0°C  isotherm.  The  front  progresses  in  a  step-like  pattern.  This  occurs 
when  the  location  of  phase  change  moves  from  one  node  to  the  next  and  is  inherent  in  the  appar¬ 
ent  heat  capacity  method  in  a  discretized  space.  Nevertheless,  the  results  of  ADIPC  show  good 
agreement  with  the  analytical  solution.  A  more  accurate  computed  solution  could  be  obtained 
by  using  a  smaller  internodal  spacing  and  smaller  time  step.  A  copy  of  ADIPC  and  output  for  this 
solution  is  included  in  Appendix  B. 


*  1  I  » _ I _ 1  -J 

0  12  3 


Time  (hrs) 


Figure  20.  Comparison  of  ADIPC  with  Neu¬ 
mann  solution. 


Comparison  of  ADIPC  with  analytical  results -two-dimensional  phase  change  verification 

Next  compare  the  results  of  ADIPC  with  an  analytical  solution  involving  phase  change  around 
a  pipe.  This  problem  is  two-dimensional  in  Cartesian  coordinates  and  thus  is  a  two-dimensional 
verification  of  ADIPC.  But  it  is  a  one-dimensional  problem  in  cylindrical  coordinates,  facilitating 
an  analytical  solution.  An  exact  analytical  solution  is  available  for  the  case  of  freezing  in  a  region, 
initially  at  a  uniform  temperature,  that  is  suddenly  subject  to  the  effects  of  a  continuous  line 
source  that  extracts  heat  at  a  rate  of  Q  per  unit  time.  The  exact  solution  of  this  heat  conduction 
problem  is  given  in  Carslaw  and  Jaeger  (1959).  The  results  are  as  follows. 

The  location  of  the  freezing  front  at  any  time  is  given  by 

R  =  2\sM[t  (84) 

where  R  is  the  radius  of  the  front  (r  -  0  is  the  location  of  the  line  source),  a,  is  the  thermal  dif- 
fusivity  of  the  frozen  zone,  t  is  time  and  X  is  obtained  from  the  following  relation: 

£  exp  ("X2)  +  E,Wa,/aJ  exp  (‘X*  °i /o*>  =  ^xiLP  <85> 


where  k  =  conductivity 
a  =  diffusivity 
-  initial  temperature 
T(  =  fusion  temperature 
/.p  =  latent  heat  per  unit  volume. 

E,  is  the  exponential  intcrgral  function 

..  .  ,  f  x  evdV 

and  the  subscripts  1  and  2  represent  the  frozen  and  unfrozen  zones,  respectively. 
The  temperatures  in  the  frozen  and  unfrozen  zones  arc  given  by 

M-Jv)  0<"K 


(86) 


27 


$ 

f' 


M 

9 


i 


> 


Ti~T{ 


(T{  -  Tf) 

Ei  (-X2  aja2) 


r  >  R  . 


(87) 


In  order  to  get  a  solution  for  phase  change  around  a  pipe  from  this,  the  pipe  is  assigned  a  con¬ 
stant  radius  rp;  the  temperature  at  rp  varies  with  time  according  to  eq  86.  A  shifted  time  is  used 
so  that  at  time  t  =  0  in  the  computer  run,  the  location  of  the  phase  change  front  is  at  r  =  r  .  For 
this  time,  the  initial  temperature  distribution  is  figured  from  eq  87. 

For  the  comparison,  the  following  values  were  used: 

k.  =  0.0072  -Cal0^  k2  =  0.0042  -C--5r  ■■ 

1  cm  s  C  z  cm  s  C 

a,  =0.014165  —  a,  =0.005556  — 

1  s  s 

Lp  =  33.012  — - 

cm2  °C 

=  4°C 

Tf  =  0°C  . 

Densities  for  all  regions  were  assumed  equal  in  the  solution. 

Equation  84  was  solved  by  substituting  a  polynomial  approximation  for  the  exponential  inte¬ 
grals  (Abramowitz  and  Stegun  1970)  and  then  solving  the  equation  for  X  by  an  iterative  scheme 
on  the  computer.  For  this  case  it  was  found  that  X  =  0.08246. 

ADIPC  used  an  internodal  distance  of  l  cm  and  a  time  step  of  60  s  on  a  30  by  30  grid.  Because 
of  the  symmetry  of  the  problem,  with  the  Cartesian  origin  at  the  center  of  the  pipe,  only  one- 
quarter  of  the  situation  was  modeled.  The  top  and  right-hand  edges  of  the  grid  were  assigned  zero 
flux  boundaries,  and  the  left-hand  and  bottom  edges  of  the  grid  were  assigned  the  semi-infinite 
condition. 


Figure  21.  Comparison  of  ADIPC  with  radial  analytic 
freezing  solution. 
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The  location  of  the  freezing  front  over  time  was  found  in  ADIPC  by  interpolating  to  find  the 
zero  degree  isotherm;  this  is  plotted  against  the  analytical  solution  in  Figure  21 .  Excellent  agree¬ 
ment  is  found. 


USER  INSTRUCTION  FOR  ADI 

The  first  step  in  using  ADI  to  solve  a  problem  is  to  define  the  boundaries  of  the  problem  and 
to  identify  the  different  materials  in  the  problem.  It  is  important  to  know  the  dimensions  of  ob¬ 
jects  to  be  modeled  as  accurately  as  possible.  Next,  identify  the  boundary  conditions  in  the  prob¬ 
lem.  Now  draw  the  grid,  determining  a  nodal  spacing  that  will  accurately  represent  the  materials 
in  the  problem.  Indicate  the  boundary  conditions  on  the  grid,  then  set  these  conditions  in  the 
computer  program  by  assigning  the  appropriate  values  of  RAY  (I,  J,  2)  to  each  node  I,  J  in  the 
grid;  the  values  are  listed  in  the  comment  statements  at  the  beginning  of  subroutine  ADDATA. 

Material  properties  such  as  thermal  conductivity,  density,  etc.,  must  be  known  in  consistent 
units.  There  are  no  dimensional  constants  inherent  in  the  program,  so  any  system  may  be  used. 

For  the  problems  presented  in  this  report,  the  following  units  were  used:  thermal  conductivity 
(W/m  K),  density  (kg/m3),  specific  heat  (W  hr/kg  K),  temperature  (°C),  time  (hr)  and  distance  (m). 
The  units  of  the  convection  coefficient  would  be  W/m2K.  Values  for  the  conductivity,  specific 
heat,  density  and  convection  coefficient  must  be  specified  by  the  user  as  indicated  in  the  main 
program  (near  the  start  of  the  2001  and  2002  loops).  These  values  may  be  programmed  to  change 
with  temperature  or  time  and  could  conceivably  be  different  for  each  node  in  the  grid.  At  the 
time  of  this  report’s  publication,  the  conductivity  of  a  semi-infinite  node  in  the  program  is  assumed 
equal  to  that  of  the  adjacent  regular  interior  node. 

Now  edit  subroutine  ADDATA  following  the  directions  in  the  comment  statements  in  the  sub¬ 
routine  to  initialize  the  variables  and  arrays.  Note  that  the  variables  and  arrays  are  defined  at  the 
beginning  of  ADDATA. 

When  editing  subroutine  ADDATA,  the  user  will  encounter  a  variable  named  ITRT,  which  is 
defined  as  the  number  of  time  steps  before  the  results  are  printed.  Setting  this  value  as  five,  for 
example,  will  set  a  counter  in  the  main  program  so  that  at  every  fifth  complete  time  step  the  main 
program  will  print  the  temperature  distribution  (RAY  [I,  J,  1  ] )  into  file  ADITMP.  Similarly,  var¬ 
iable  ITPC  will  set  a  counter  to  call  subroutine  ISOTHM  to  locate  the  user-specified  isotherms  in 
the  current  temperature  distribution.  The  coordinates  of  these  isotherms  are  printed  into  file 
POINTS. 

When  ADI  is  run,  it  will  interactively  ask  whether  or  not  the  user  wishes  to  run  subroutine 
ADDATA.  The  first  time  the  program  is  run,  ADDATA  must  be  run.  The  variables  will  be  put 
into  a  formatted  data  file,  AD  ID  AT.  If,  after  the  program  is  run,  a  change  is  made  in  the  main 
program  but  not  to  ADDATA,  the  subroutine  need  not  be  called  again;  the  main  program  will 
read  the  data  from  ADIDAT. 

Once  the  program  has  been  run,  the  user  should  examine  output  file  ADIOUT  closely  to  be 
sure  that  the  initial  conditions  and  boundary  conditions  are  those  intended,  i.e.,  that  no  mistake 
was  made  when  editing  ADIDAT  or  changing  values  in  the  main  program. 

USER  INSTRUCTION  FOR  ADIPC 

The  use  of  ADIPC  is  the  same  as  ADI  with  the  exception  of  the  specification  of  several  phase 
change  variables.  When  editing  ADDATA,  the  user  must  specify  the  value  of  the  latent  heat  of 
fusion,  the  conductivity  and  density  of  the  phase  change  region,  and  the  temperature  range  over 
which  phase  change  will  occur.  If  the  units  indicated  in  the  instructions  for  ADI  are  used,  the 
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latent  heat  will  have  the  units  W  hr/kg,  as  does  the  specific  heat.  The  value  of  the  apparent  specific 
heat  (discussed  in  the  Phase  Change  section),  known  as  CPPC  in  the  program,  should  be  calculated  by 
the  user  as  follows: 


CPPC=  lA  (CpB 


where  CPb  =  specific  heat  of  the  frozen  material 

CpA  =  specific  heat  of  the  unfrozen  material 
Hl  =  latent  heat  of  fusion  for  the  material 

AT  =  temperature  range  around  the  fusion  temperature  where  phase  change  occurs  from 


T<- 


A  T  AT 

T  to7>+T 


AT  in  the  verifications  was  1 .0°C. 

As  in  ADI,  the  values  for  density,  specific  heat  and  conductivity  must  be  specified  as  indicated 
in  the  comments  near  the  start  of  the  main  program,  ADIPC. 


CONCLUSIONS 

Two  two-dimensional  finite  difference  computer  programs  have  been  developed  to  model  time- 
varying  heat  conduction.  Results  of  test  runs  of  the  programs  show  excellent  agreement  with  an¬ 
alytical  and  experimental  results.  The  programs  are  easily  set  up  to  model  new  problems,  and  have 
the  capability  to  solve  a  wide  variety  of  heat  conduction  problems. 
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APPENDIX  A:  PROGRAM  INFSUM  AND  SAMPLE  INPUT 
AND  OUTPUT  FOR  PROGRAM  ADI 


INFSUM,  a  short  program  to  do  infinite  sum  calculation 
for  rectangle  problem 


c 


INFSUM 

TrilS  CALCULATE  S  Til!  INFINITE  '>JM  PRODUCT 
IMPLICIT  DOUBLE  P  K  F  •"  1 3  I  0  N  t  A-H*0-2> 
IMPLICIT  INTCGtR«<t<  l-*ii) 

UImLNSION  T  E  /.  P  (  N  J  »  d  -  ) 

CALL  C3NTRL  <2  ,  •  IxlFOJT  ♦  .5) 

*,-•  —  •  i 


FOR 


THE  RECTANGLE. 


C 


220 

221 


23G 

231 

210 

20G 


211 

212 

103 


A  =  .  -  fl  2  5 1 
LL-i 

ri  - 

T  u -ID  * 

T  1  =  2  C  •  0 

1  I-^l'tlH02fc.:4 
P  I  2=  1.1,093  J  44  C4 
[■  0“  •  1 
JCLT=1 
I  »  =  i.  1 
I  r  x  2  \ 

i.  0  ICC  1  T  =  b  ,  M  «  j 

T  1H-I  T  OEL  T 
JO  CFO  1  =  1,  IT 
JO  210  J  =  1  ,  I  X 
r=L.i»<  i - 1 .  > 

X=JL •< J-l, ) 

Cl  =  -A»Pf2*Tl;T/(>».*EL*FLl 
C2=PI*X/(2.»EL) 

C3  =  -A*P  [  ?»T  J 

C4  =  F'I*Y/(2.«LO 

FIOUKr.  SUM 

Ns  1 

3  L‘  •  1 1  =  ^  . 

JO  2  2  0  I  \  N  =  1  ,  L 
I  i  -  I  \l  N  - 1 

CIA  =  C1*<2.*1N*1. >*<•'. *IL»1.> 

C  1  A  A  =  f)  t  A  F  (CIA) 

C?A=C2« <2. *  I N ♦  1  •  > 

FH 0 ( =< (N/(2 .• IN  +  l  .  )  )»C 1 AA  •  CO  ■  1C2A  )  ) 
SU-'1=SJM1  +  FKST 
N  =  -  i  •  N 

IF<C1AA.LE..0CC  1J  G ..  TO  221 

CONTINUE 

CONTINUE 


FIGUrtC  S.UM2 

N  =  1 

$  J'-*2  =  3. 

no  23C  I  N  N  =  1 , 5  0 
I  \  =  I NN-  1 

CJA=C3»(2»*IN*1.)*I'.»JN*1.> 

C4A  =  C4*  <2. *  I N  ♦  1  »  ) 

C3AA=UCXP(C3A) 

SCNt>s(N/(2#*IK*l«))*C3‘A*CC0S(C4A) 

SUM?=SUM2*SCN0 

N=-1*N 

IF(  C36A  .LL  .  .3  00  1  )  C,„  TO  23  i 

CONTINUE 

CONTINUE 


T  t  M  P  C I ,J)=T  OMTI-TC  >  * ( 1  fc  •  /  P I  2  >  *  s  U  M 1 *SUM2 

CONTINUE 

CONTINUE 

GRITE(b,211>  TIM 

FORMAT  (  /  , IX,* TEMPERATURES  AFTER»,Ff,.2,* 
WRITE(b,212)  («TEMP(I,J|,J=1,IX»,I=1,IV) 
FORMAT  (1X«11F6.2) 

CONTINUE 

CALL  C3NTRL(4»*INF0JT*,5) 

CALL  EXIT 
END 


HOURS*  ) 
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rr 

-  -  -  -_v 

O  • 

« . 

*  *  * 

? 

Sample  output  from  program  INFSUM 

>Jk  * 

TEMPERATURES 

AFTER 

3.0  0 

HOURS 

20. CO 

20.00 

20.30 

23.00 

20  .CO 

20.  r  3 

20  .23 

21.17 

24.14 

30.56 

40.00 

►  .*■ 

20. 0C 

2C.0C 

20.03 

23.00 

20.00 

20.03 

20.23 

21.17 

24. :  4 

30.56 

40.00 

f  *  * 
*>  * 

20.00 

20.00 

20.00 

20.00 

20.00 

20.03 

20.23 

21.17 

24.14 

30.56 

40.00 

wN 

20.00 

20.00 

20.03 

2  0.00 

2  0.00 

20.03 

20.23 

21  .17 

24.14 

30.56 

40  .00 

n 

20.00 

2  0.00 

20.00 

2  L  .  0  0 

20.00 

23.03 

20.23 

21.17 

24.14 

30.5b 

40.00 

20.00 

23.00 

20.00 

2  0 . 0  0 

20.00 

20  .  03 

20.23 

21.17 

24.14 

30.56 

40.00 

F7 

20.00 

20. 0C 

20.3  3 

2  0.03 

20.00 

20.03 

20.23 

21.17 

24.14 

30.56 

40.00 

fv 

20.00 

20.0  3 

20.33 

2  0.0  n 

20.00 

20.02 

20  .23 

21.17 

24.14 

30  .56 

40  .00 

tv 

20.00 

2  0.00 

2  0.00 

23.00 

20.00 

20.03 

20  .23 

01.17 

24.14 

3  0.56 

40.00 

r  " 

20.00 

2  0.00 

20.00 

2C.00 

20.00 

20.03 

20.23 

21.17 

24.14 

30.56 

40.00 

F-\ 

20.00 

23.00 

20.00 

20.30 

20.03 

20.03 

20.23 

21.17 

24.14 

30.56 

40.00 

r**- 

20.00 

2  3 . 0  0 

20.00 

2  0.03 

20.00 

2  0.03 

20  .23 

21.17 

24.14 

30.56 

40  .00 

r«"i 

.'0.0  0 

2  3 . 0  J 

2U.03 

21.00 

20. CO 

20.03 

20.23 

21.17 

24.14 

30.56 

40.00 

[V, 

2  0.00 

2  0.03 

20.03 

2  u  .  0  C 

2  0.00 

20.03 

20.23 

21.17 

24.14 

30.56 

40.00 

H 

20.00 

2  3.00 

20.00 

2  u  .  0  0 

20.01 

20.04 

20.23 

21  .17 

24.14 

30.56 

40.00 

20.03 

2  0.03 

20. C3 

2  3.33 

20.04 

20.06 

20  .26 

21  .20 

24.16 

30.57 

40  .00 

V 

20  .23 

v  J  .23 

20.23 

20.23 

2  0.23 

2  0.26 

20.46 

21  .38 

24.32 

30.67 

40.00 

-  T* 

21.17 

21.17 

21.17 

21.17 

21.17 

21.20 

21  .38 

22.26 

25.06 

31.11 

40.00 

v'' 

2<*  .  1 4 

24.14 

24.14 

24.14 

24.14 

24.16 

24 .32 

25.06 

27.42 

32.51 

40.00 

a© 

30  .do 

3  J  .56 

3  0.56 

3  0.56 

3  J  .56 

30.57 

30  .67 

31  .11 

32.51 

35.54 

40.00 

Kj 

40.00 

4  J  .  0  0 

40.00 

4  J  .  3  3 

40.00 

40.00 

40.0  0 

40. oc 

40.00 

40.00 

40.00 

m 

TEMPER  ATOKtS 

AFTER 

K.  00 

HOURS 

23  .JO 

23.03 

20.01 

23.04 

20.15 

20.51 

21.4b 

23.61 

27.44 

33.11 

40.00 

2  0  .  0  0 

2  0.03 

20.01 

2  C  .  0  4 

20.15 

20.51 

21.48 

23.61 

27.44 

33.11 

40.00 

20.00 

2G.00 

20. 31 

2  0.04 

2  0.15 

20.51 

21.48 

23.61 

27.44 

33.11 

40.00 

JT 

20.00 

2  3.00 

23.01 

2  *  •  G  4 

20.15 

20.51 

21.48 

23.61 

27.44 

33.11 

40.00 

20.03 

2  3.00 

20.01 

2  3.04 

2C.15 

2U.51 

21.48 

23.61 

27.44 

33.11 

40.00 

23.00 

2  0.00 

20.31 

2  0.04 

20.15 

20.51 

21.48 

23.61 

27.44 

33.11 

40.00 

20.0  0 

j  0 . 0  0 

2  0  .  ul 

2  3.04 

2  0.16 

2  0.51 

21.48 

23.61 

27.44 

33.11 

40.00 

20.30 

2  0.00 

20.01 

2  3.34 

20.15 

20.51 

21.48 

23.61 

27.44 

33.11 

40.00 

23  .  j  0 

2  0.03 

2  3  .  J 1 

2  3.04 

20.15 

2  0.51 

21.48 

23.61 

27.44 

33.11 

40.00 

20.00 

2  j  .  Q  3 

20.01 

2  3.04 

20.15 

23.51 

21.48 

23.61 

27.44 

33.11 

40.00 

20 .00 

2  3.03 

20.31 

20. 34 

20.15 

2  0.51 

21.48 

23.61 

27.44 

33.11 

40.00 

2C.0Q 

2  0.00 

20.01 

20.04 

20 . 15 

20.51 

21.49 

23.61 

27.44 

33.11 

40.00 

23.01 

2  3.01 

20.31 

2  2.04 

20.16 

2  0.52 

21.49 

23.62 

27.45 

33.11 

40.00 

2  0 . 3  4 

2  0.04 

2  0  .  u  4 

2  t  .  0  7 

20.18 

2'  3 . 5  5 

21  .52 

23.64 

27.46 

33.12 

40.00 

m 

20. Id 

2  U  •  1  d 

20.16 

2  C  .  1  ft 

2G.30 

j  0 . 6  6 

21.62 

23.73 

27.53 

33.16 

40.00 

PLn 

2  0.61 

20.51 

2  0.52 

2  0.55 

20.66 

21.01 

21.96 

24.03 

27.76 

33.28 

40.00 

21  ,4rt 

21.44 

21.49 

21.52 

21.62 

21.96 

22.86 

24.83 

28.37 

33.62 

40  .00 

1 1 

23. M 

2  3.61 

2  3.62 

2  3.64 

23.73 

8  4. 03 

24  .63 

26.57 

29.71 

34.35 

40  .00 

3  V 

27.44 

27.44 

2  7.45 

7.46 

27.53 

27.  76 

28.37 

29.71 

32.11 

35.67 

40.00 

SNc 

33.11 

33.11 

33.11 

33.  12 

33.16 

33.2b 

33.62 

34.35 

35.67 

37.62 

40.00 

& 

40.00 

4  0.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

TEMPERATURES 

AFTtK 

1  5 . 0  0 

HOURS 

*  ,j 

20.01 

20.02 

20.07 

2  0.21 

20.58 

21.37 

22.90 

25.49 

29.32 

34.31 

40.00 

jjt 

20.01 

20.02 

20.07 

20.21 

20.58 

21.37 

22.90 

25.49 

29.32 

34.31 

40.00 

is 

20.01 

20.02 

20.07 

20.21 

20.58 

21.37 

22.90 

25.49 

29.32 

34.31 

40.00 

20.01 

20.02 

20. 07 

20.21 

20.58 

21.37 

22.90 

25.49 

29.32 

34.31 

40.00 

20.01 

20.02 

23.07 

20.21 

20.58 

21.37 

22.90 

25.49 

29.32 

34.31 

40.00 

tv. 

20.01 

20.02 

20.07 

23.21 

20.58 

21.37 

22.90 

25.49 

29.32 

34.31 

40  .00 

w 

20.01 

20.02 

20.07 

20.21 

20.58 

21.37 

22.90 

25.49 

29.32 

34.31 

40.00 

55 

20.01 

20.02 

20.07 

20.21 

20.58 

21.37 

22.90 

25.49 

29.32 

34.31 

40 .00 

u 

20.01 

20.02 

20.07 

2  0.22 

20.58 

21.37 

22.90 

25.49 

29.32 

34.31 

40.00 

20.01 

20.02 

2  0.07 

20.22 

20.56 

21.37 

22.90 

25.49 

29.32 

34.31 

4C.00 

20.02 

20.03 

2  r  .03 

20.22 

20.56 

21.37 

22.90 

25.49 

29.32 

34.31 

40.00 

s 

20.03 

2  0.04 

2  •  j9 

20.24 

20.60 

21.39 

22.92 

25.50 

29.33 

34.32 

40  .00 

*m* 

20.08 

20.09 

2  j  .  14 

20.29 

20.64 

21.43 

22.96 

25.54 

29.36 

34.33 

40.00 

»■  , 

20.23 

20.24 

20.29 

20 . 43 

20.7b 

21.57 

23.08 

25.64 

29.44 

34.37 

40 .00 

20.09 

20.60 

20. 64 

20.78 

21.13 

21.90 

23.39 

25.90 

29.63 

34.47 

40.00 

W 

21.38 

21.39 

21.44 

21.67 

21.90 

22.64 

24  .07 

26.48 

30.05 

34.70 

40.00 

SC 

22.91 

22.92 

22.46 

23.08 

23.39 

24.07 

25.38 

27.59 

30.87 

35.14 

40.00 

w 

25.49 

25.50 

25.54 

25.64 

25.90 

26.4b 

27.59 

24.47 

32.25 

35.87 

40.00 

29.33 

29.33 

29.36 

29.44 

29.63 

30.05 

30  .87 

32.25 

34.30 

36.96 

40.00 

I  34.31 

34.32 

34.33 

34.37 

34.47 

34.70 

35.14 

35.87 

36.96 

38.38 

40 .00 

?! 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.0  0 

40.00 

40.00 

40.00 

TEMPERATURES 

AFTf  R 

23.00 

HOURS 

20.06 

20.10 

20.23 

20.54 

21.17 

22.29 

24.14 

26.87 

30.56 

35.05 

40  .00 

20  .06 

20.10 

20.23 

20.54 

21.17 

22.29 

24.14 

26.87 

30.56 

35.05 

40.00 

20.06 

20.10 

20.23 

20.54 

21.17 

22.29 

24.14 

26.5' 

30  „56 

35.05 

40.00 

!«  • 

20.06 

2C.10 

20.23 

20.54 

21.17 

22.29 

24.14 

26.8. 

30. 56 

35.05 

40.00 

20.06 

20.10 

20.23 

20.54 

21.17 

22. 29 

24.14 

26.87 

30.56 

35.05 

40.00 

20.06 

20.10 

20.23 

20.54 

21.17 

22.29 

24.14 

26.87 

30.6 

35.05 

40.00 

[  20.06 

20.10 

20.23 

20.54 

21.17 

22.29 

24.14 

26.88 

30.56 

35.05 

40.00 

4  4 

1  20.06 

20.10 

20.24 

20.54 

21.17 

22.29 

24.14 

26.88 

30.53 

35.05 

40.00 

> 

20.07 

20.10 

20.24 

20.55 

21.17 

22.29 

24.14 

26.88 

30.56 

35.05 

40.00 

Cy 

20.07 

20.11 

20.24 

20.55 

21.18 

22.30 

24.14 

26.88 

30.56 

35.05 

40.00 

LN 

20.10 

20.13 

20.27 

20.58 

21.20 

22.32 

24.16 

26.90 

30.57 

35.05 

40  .00 

rV 

20.15 

20.19 

20.32 

20.63 

21.25 

22.37 

24.21 

26.93 

30.60 

35.07 

40.00 

K 

20.29 

20.33 

20.46 

20.77 

21.38 

22.50 

24.32 

27.03 

30.67 

35.10 

40  .00 

3 

20.61 

20.64 

20.77 

21.07 

21 .68 

22.77 

24.57 

27.23 

30.81 

35.18 

40.00 

21.23 

21.26 

21.39 

21.68 

22.26 

23.32 

25.06 

27.64 

31.11 

35.33 

40 .00 

N 

22.35 

22.38 

22.50 

22.77 

23.32 

24.32 

25.95 

28.38 

31.64 

35.61 

40.00 

24.19 

24.22 

24.32 

24.57 

25.06 

25.95 

27.42 

29.59 

32.51 

36.07 

40.00 

'  26.92 

26.94 

27.03 

27.23 

27.64 

28.38 

29.59 

31.39 

33.80 

36.75 

40.00 

30.59 

30.61 

30.67 

30.82 

31.11 

31.64 

32.51 

33.80 

35.54 

37.66 

40  .00 

<£■ 

35.06 

35.07 

35.10 

35.18 

35.33 

35.61 

36.07 

36.75 

37.66 

38.77 

40.00 

I 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

32 

40.00 

40.00 

40.00 

40.00 

40.00 

1 

-V  _ 

f'-.‘  -  ~ 

|,N  .**.*•*  .  V  %  ,\V 

y-y-y- 

«  '  ^9 

.  ■  -  - 

* 

4*k4 

jvlA-'  vA 

m.  *\  *  .  O 

i  • 

* .  < 

Data  printed  by  program  ADI  for  rectangle  problem  (file  ADIOUT) 


DATA  FOR  THIS  RUN  DF  ADKj: 

OS  CELT  01 

0.10U.C  o.CCOOQ  SO. 00000 

A  X  Y  NISO  ITKT  IMAX 

1  XI  01  S  1  10 

T ISO(H>  »H=l,MISO: 

40.00 

30.00 

32.00 

28.00 

24.00 


<  RNR < 1 )  : 


KRNR<2>=  1  K*’lR(3>  =  1  KkNR{4>= 


.XT  < J) 

F  LXy< J) 

TMPT  ( J  ) 

T1PB ( J) 

0 . 0  C  3  3 

o.oroo 

?  0 . 0  0  0  c 

40.000C 

0.0000 

o .  c  o  a  o 

20.0000 

40.0000 

o.ccno 

u .  o  :•  c  o 

20.0000 

40.0000 

0.0000 

C .0  100 

20.0003 

40.0000 

0.0030 

0.0003 

20.0000 

40.0000 

P.OCC j 

0 . C  10  0 

20.0000 

40.0000 

J  .  0  0  0  0 

3.000 

20.0000 

40.0000 

u.gcoo 

0  .  0  ..  C  0 

2  0 . 3  0  0  C 

40.0000 

3.0000 

0  .  u  J  C  0 

20.0000 

40.0000 

0.0000 

C  .  0  1 0  0 

20.0000 

40.0000 

■j  .  0  0  j  0 

O.C "03 

20.0000 

40.0000 

C.OOC  3 

J  .  1  .<  0  0 

20.0000 

0 .030  J 

C  .  r  :  5  0 

20.3000 

0  .  000  0 

o  .  Cl 0  0 

2C.Q030 

0 . 0  C  0  0 

0.0300 

20.00  OC 

0.000  0 

0 . 0  0  0 

2  0.0000 

C  .000  0 

0.000 

2U.C0DC 

0  «  C  0  0  J 

0  .  j  J  0  0 

20.0000 

J  .  0  0  J  0 

o.oroo 

2  0  .  C  0  0  C 

0.0QJ0 

0.0:00 

20.0000 

o .  c  o  o : 

0.0300 

20.0000 

0.0000 

0.0003 

20.0000 

0  .  0  C  0  3 

0.0100 

20.0000 

0 . 0  C  0  0 

0.0000 

20.0000 

0  .  C  ij  0  0 

o.oooo 

20.0000 

3 . 0  0  u  0 

0.0000 

20.0000 

o.OUj 

O.C  -  0  u 

23.0000 

o.oooo 

0.0000 

20.0000 

u  .  o  C  0  0 

0 . 0  '  0  0 

20. COCO 

0.000  0 

0.0000 

20.0000 

0.  000  0 

O.C coo 

2  0  .  0  0  0  C 

0.0030 

0.0000 

20.0000 

4  A  Y (  I,J 

,*)  NPQAL  MATERIAL 

TYPE 

1.00 

1.00 

l.CC 

1 .00 

1.00 

1  .  oc 

1.00 

1  ,G0 

1.0  0 

l.oa 

1.00 

1.00 

1.00 

1.00 

1  .CO 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.C0 

1.0c 

1.00 

1.0c 

1.00 

1.00 

1.3  J 

l.oo 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.  CO 

1.00 

1.00 

1.00 

1.00 

1.03 

1.30 

1.00 

1.00 

1.00 

1.00 

1  .00 

1.00 

1.00 

1.00 

1.30 

1.00 

1.00 

1.00 

1.00 

1.00 

1.0  0 

1.00 

1.00 

1.0  0 

1.00 

1.03 

1 .  C  3 

1.00 

1.00 

1.00 

1.0c 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.05 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.C0 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.03 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1  .03 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.  00 

1.00 

1.00 

1.00 

1.00 

R  A  Y  <  I  »  U  «  2  )  NODAL  LOCATION 

TYPE 

00 

8 

00 

5 

00 

s 

00 

5 

00 

5 

00 

5 

00 

5 

00 

5 

00 

s 

00 

5 

40.0000 
40.0000 
40.0000 
40.0000 
40.0000 
40.0000 
40.0000 
40.0000 
40.0000 
4  0.0  00  0 
40.0000 
40.0000 
4  0.000  0 
40.0000 
40.0000 
40.0000 
40.0000 
40.0000 
40.0000 
40.0000 
40.0000 


00 

8 

0 

0  8 

.00 

8 

00 

5 

0 

0  5 

.00 

5 

00 

5 

0 

0  5 

00 

5 

00 

5 

0 

0  5 

00 

5 

00 

s 

0 

0  S 

00 

5 

00 

5 

0 

0  5 

00 

5 

00 

5 

0 

0  5 

00 

5 

00 

5 

0 

0  5 

00 

5 

00 

5 

0 

0  5 

00 

5 

00 

5 

0 

0  5 

00 

5 

V 


3 

.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

7.00 

8 

.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

7.00 

8 

.00 

5.CG 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

7.00 

8 

.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

7.00 

8 

.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

7.00 

8 

.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

7.00 

9 

.00 

5.00 

5.00 

5.0  0 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

7.00 

8 

.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

7.00 

3 

.00 

5.00 

3.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

7.00 

9 

.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

7.00 

2 

.00 

7.00 

7.00 

7.00 

7.00 

7.00 

7.00 

7.00 

7.00 

7.00 

3.00 

RAYtliJ. 

1)  TEMPERATURES  AT 

THE  START: 

20 

.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

20.03 

20.00 

20.00 

20.00 

20.00 

20.00 

20.  CQ 

20.00 

20.00 

AO. 00 

20 

.00 

20.00 

20. 0C 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

20.00 

20. 00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

•  u  0 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20. CO 

AO. 00 

20 

.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

20.00 

20.00 

20.00 

20. OC 

20.00 

20.00 

20.00 

20.00 

20. CO 

AO. 00 

20 

.00 

20. CC 

20.  ac 

20. 0C 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

20.00 

20.0  0 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

23 

.00 

20.00 

20. 3C 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

20.00 

2  0.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

2  0.00 

20.00 

20 .00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

2  0.00 

2C.0J 

20.00 

20.30 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

2  0.00 

20.00 

20.0  0 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

2  0.00 

2  C  .  0  0 

20.00 

20. CO 

20.00 

20.00 

20.0  0 

20.00 

20.00 

AO. 00 

20 

.00 

20.00 

20.00 

20. CO 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

20 

.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

20.00 

AO. 00 

A0 

.00 

A0. CO 

A0. 00 

AO. 00 

AO  .00 

AO. 00 

AC  .00 

AO. 00 

AO. 00 

AO. 00 

AO. 00 

TEMPERATURE  S  WILL  BE  PRINTED  EVERY  1  TIME  STEPS. 

BOUNDARY  CONDITIONS: 

TOP  DTjNDARY  CONSTANT  HEAT  FLUX 

TOP  LEFT  CORNER  CONST  FLUX  ON  30TH  SIDES 

LEFT  jOUNJARY  CONSTANT  HEAT  FLUX 

TOTTOM  LIFT  CORNER  CONSTANT  TEMPERATURE 
BOTTOM  BOJNOARY  CONSTANT  TEMPERATURE 
BOTTOM  RIGHT  CORNERCONSTANT  TEMPERATURE 
RIGHT  bOUNDARY  CONSTANT  TEMPERATURE 

TOP  RIGHT  CORNER  CONSTANT  TEMPERATURE 

NUMBER  OF  TIME  STEPS=  10 


Sample  output  from  program  ADI  for  rectangle  problem  (file  ADITMP) 


R  AY  U 

.  J  . 

1) 

TEMPERATURES 

AT 

THE 

start: 

20. 

00 

20. 

00 

23  . 

00 

20.0  0 

23 

.00 

23 

.00 

20.0  0 

20. 

00 

20. 

00 

20.00 

20 

.00 

23 

.00 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20 

.00 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20, 

.0  0 

20.00 

20. 

00 

20. 

30 

20. 

00 

20.00 

20 

.00 

20 

.0  0 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.03 

20 

.0  0 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

2  C 

.00 

20 

.00 

20.00 

20. 

00 

20. 

00 

20  . 

00 

20.00 

20 

.00 

20 

.0  0 

20.00 

20. 

00 

20. 

00 

20  . 

00 

20.00 

20 

.00 

20 

.00 

20.00 

20  . 

0  0 

20. 

00 

20. 

00 

20.00 

2  U 

.  0  3 

20 

.0  0 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.  0  G 

20 

.0  0 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20 

.00 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20 

.  0  0 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20 

.00 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20 

.0  0 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20 

.00 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20 

.00 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20 

.0  0 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20 

.0  0 

20.00 

20. 

00 

20. 

00 

20. 

00 

20.00 

20 

.00 

20 

.00 

20.00 

20. 

0  0 

AO. 

OQ 

AO. 

00 

AO  .00 

AO 

.00 

AO 

.00 

AO  .00 

AO. 

00 

TEMPERATURES 

AFTER  1 

TIME  : 

STEPS  ( 

5.00  HRS): 

3A 

20. 

00 

20. 

00 

20.00 

20 

.01 

20 

.03 

20.10 

2  U  • 

20. 

00 

20. 

00 

20.  QQ 

20 

.01 

20 

.03 

20.10 

20. 

3A 

20. 

00 

20. 

00 

20.00 

20 

.01 

20 

.03 

20.10 

20. 

3  A 

20. 

00 

20. 

00 

20.00 

20 

.01 

20 

.03 

20.10 

20. 

3A 

20. 

00 

20. 

00 

20.00 

20 

.01 

20 

.03 

20.10 

20. 

3A 

20. 

00 

20. 

00 

20.00 

20 

.01 

20 

.03 

20.10 

20. 

3A 

20. 

00 

20. 

00 

20.00 

20 

.01 

20 

.03 

20  .10 

20. 

3  A 

20. 

00 

20. 

00 

20.00 

20 

.01 

20 

.03 

20.10 

20. 

3A 

20. 

.00 

20. 

00 

20.00 

20 

.01 

20 

.03 

20.10 

20. 

3A 

20. 

00 

20. 

00 

20.00 

20 

.01 

20 

.03 

20.10 

20. 

3A 

20. 

.00 

20. 

.00 

20  .00 

20 

.01 

20 

.03 

20.10 

20. 

.3A 

20. 

.00 

20. 

.00 

20.00 

20 

.01 

20 

.03 

20.10 

20. 

,  3  A 

20. 

.00 

20. 

.00 

20.01 

20 

.01 

20 

.03 

20.11 

20. 

,  3A 

34 


20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

2  0. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20, 

00 

20. 

00 

20. 

00 

AO 

.00 

20, 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

00 

20. 

00 

20. 

00 

AO 

.00 

20. 

.00 

20. 

00 

20. 

00 

AO 

.00 

20. 

,00 

20. 

00 

20. 

00 

AO 

.00 

20 

,00 

20. 

00 

20. 

00 

AO 

.00 

20. 

.00 

20. 

00 

20. 

00 

AO 

.00 

AO 

,00 

AO 

.00 

AO. 

00 

AO 

.00 

21 

.12 

23. 

70 

32. 

16 

AO 

.00 

21 

,12 

23. 

.70 

32. 

16 

AO 

.00 

21 

.12 

23 

.70 

32. 

16 

AO 

.00 

21 

.12 

23 

.70 

32. 

16 

AO 

.00 

21 

.12 

23 

.70 

32. 

16 

AO 

.00 

21 

.12 

23 

.70 

32. 

16 

AO 

.00 

21 

.12 

23 

.70 

32. 

16 

AO 

.00 

21 

.12 

23 

.70 

32. 

16 

AO 

.00 

21 

.12 

23 

.70 

32. 

16 

AO 

.00 

21 

.12 

23 

.70 

32. 

16 

AO 

.00 

21 

.12 

23 

.70 

32. 

16 

AO 

.00 

21 

.13 

23 

.70 

32. 

16 

AO 

.00 

21 

.13 

23 

.70 

32. 

16 

AO 

.00 

20.01 

20.01 

20.01 

2  J  .  02 

20.04 

20.11 

20.35 

21.13 

23.71 

32.17 

40.00 

20.03 

20.03 

20.03 

2  0.34 

2  0.06 

20.14 

20.37 

21.15 

23.72 

32.17 

40.00 

20.10 

20.10 

20.11 

2  0.11 

23.14 

20.21 

20.44 

21.22 

23.78 

32.20 

40.00 

2  0.34 

20.34 

20.34 

2c  .56 

20.37 

20.44 

20.66 

21.45 

23.98 

32.30 

40.00 

21.12 

21.13 

21.15 

21.15 

21.16 

21.22 

21.45 

22.19 

24.61 

32.60 

40.00 

23.70 

23.70 

23.72 

2  5.71 

23.72 

23.78 

23.98 

24.61 

26.71 

33.61 

40.00 

32.16 

32.16 

32.16 

3  2 . 1  7 

32.17 

32.20 

32.39 

32.60 

33.61 

36.93 

40.00 

4  0.00 

40.00 

4  0.00 

4  0.00 

4C.0  0 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

TEMPERA? JRES  AFTER  ? 

TIME  steps  < 

10.00  HRS): 

20.01 

20.01 

2  0.03 

*2  0  .  U  7 

20.2  0 

20.55 

21.46 

23.60 

27.89 

32.98 

40.00 

20.01 

20.01 

2  0  •  u  3 

20.07 

20 .2  0 

20.55 

21.4b 

27.60 

2  7.89 

32.98 

40.00 

20.01 

20.01 

2  0  .  u  3 

2  0.07 

2  0.20 

20.55 

21.46 

23.60 

27.89 

32.98 

40.00 

20.01 

20.01 

20. 3T 

2  J  .  0  7 

20.2  0 

20.55 

21.46 

23.60 

27.89 

32.98 

40.00 

20.01 

20.01 

2  3  .  u  3 

23.0  7 

23.20 

20.55 

21.46 

23.69 

27.89 

32.98 

40.00 

20.01 

20.31 

2  0.33 

2  3.37 

2  0.20 

20.55 

21.46 

23.60 

27.89 

32.98 

40.00 

20.01 

20.01 

2  0.33 

20.37 

20.20 

20.55 

21.46 

23.69 

27.89 

32.98 

40.00 

20.01 

2  0.01 

20.0  5 

23.0  7 

20.20 

20.65 

21.45 

23.60 

27.89 

32.98 

40.00 

20.01 

20.01 

20  .3  5 

2c  .0  7 

20 .20 

20.56 

21.46 

23.60 

27.89 

32.98 

40.00 

20.01 

20.01 

20.0  3 

2  5  .  C  7 

20.2  0 

2  0.36 

21.46 

23.60 

27.89 

32.98 

40.00 

23.01 

20.01 

2  C  .  0  5 

2 1  .  C  7 

2  0.21 

20.56 

21.46 

23.60 

27.90 

32.98 

40.00 

20. til 

20.02 

20.03 

2  0.08 

2  0.21 

20.56 

21.47 

23.61 

27.90 

32.98 

40.00 

20.03 

20.03 

20  .f  3 

L  J  .  1  J 

2  J.23 

20 .58 

21.48 

23.62 

27.91 

32.99 

40.00 

2  0.08 

2  0.08 

2  0.10 

8  J  .  1  4 

20.27 

20.62 

21.53 

23.66 

27.94 

33.01 

40.00 

2C.  21 

20.21 

20 .23 

£  ...27 

23. 4C 

80.75 

21.65 

23.77 

28.02 

33.05 

40.00 

2  0.56 

2  0.56 

2  0.5  c 

.62 

20.75 

21.09 

21.97 

24.05 

28.23 

33.13 

40.00 

21.46 

21.47 

21  .46 

81 .53 

21.65 

21.97 

22.81 

24.80 

28.78 

33.49 

40.00 

2  3.40 

2  3  .  v  1 

2  3.6  V 

25.  =  6 

2  5.77 

24.05 

24.89 

26.55 

30.07 

34.24 

40.00 

27.40 

27. 

?7.il 

2  7.54 

26.0  2 

28.23 

2  8.78 

■30.07 

32.67 

35.75 

40.00 

32.96 

32.98 

5  2.3  .» 

5  5.31 

33.09 

33.18 

33.49 

34.24 

35.75 

37.54 

40.00 

4  0.00 

4  0.00 

4  0.02 

4.  .03 

4  0.00 

4  0.00 

40.00 

40.09 

49.00 

40.00 

40.00 

TEMPERATURE  S  AFTER  3 

TIME  S 

TEPS  ( 

15. 0C  1 

HRS  ) : 

20.03 

£.  .  0  4 

2  0.1 ; 

2  0  .26 

2  0.61 

21.39 

22.94 

25.60 

29.30 

34.45 

40.00 

20  .0  J 

2  0.04 

2  u  .  1  . 

2  o  .  2  u 

2  3.61 

21.39 

22.94 

25.69 

29.30 

34.45 

40.00 

20.03 

2  0.34 

20.1. 

8  .  2  6, 

20.61 

21.39 

22.94 

25.69 

29.30 

34.45 

40.00 

2  C  .  0  3 

c  0  .  O’  4 

2  0.12 

8  o  .  2  5 

20.6  1 

21.39 

22.54 

25.60 

29.30 

34.45 

40.00 

2  0.03 

2  C  .  -  4 

2  0.1 

0  0.2  ‘.i 

20.61 

81.79 

22.^4 

25.60 

29.30 

34.45 

40.00 

2  0  .  C  3 

20.04 

£>3.1 

8  ■  .25 

23.61 

21.39 

22.94 

25.60 

29.30 

34.45 

40.00 

20.03 

2  J  •  j  4 

2  0.10 

2  c  .26 

2  0.61 

21.39 

22.94 

25.60 

29.30 

34.45 

40.00 

2  0.03 

20.04 

2  0.1c 

8  2.23 

20.61 

21.39 

22.94 

25.60 

29.30 

34.45 

40.00 

2  0.03 

20.0-1 

2  0.1  , 

8  .26 

20.61 

21.39 

22.94 

25.60 

29.30 

34.45 

40.00 

2  j.03 

20 . 03 

2  0.11 

2  ...  2 6 

20.61 

21.39 

2  2.94 

25.60 

29.30 

34.45 

40.00 

20.04 

8  0.  J  ft 

2  0.12 

20.27 

20.62 

21.40 

22.95 

25.61 

29.31 

34.45 

40.00 

20.07 

20  .  Jft 

20.14 

2  0 . 2  9 

2  0.65 

21.42 

22.97 

25.63 

29.32 

34.46 

40.00 

20.13 

20.14 

20.2  1 

?  0.3ft 

20.70 

8  1.48 

23.92 

25.67 

29.35 

34.48 

40.00 

2  0 . 2  h 

20.3  0 

2  C  .  3  3 

2  v  •  5  0 

23.85 

21.62 

23.15 

25.76 

29.44 

34.52 

40.00 

20.64 

2  0.65 

20.71 

2  0.95 

21.20 

21.95 

23.45 

26.04 

29.63 

34.62 

40.00 

21.41 

21.43 

21.45 

£1.62 

21.95 

2  2.68 

24.12 

26.60 

30.04 

34.83 

40.00 

22.9ft 

22 . 9  7 

23.0  2 

2  5.15 

2  3.4'i 

24.12 

25.44 

27.71 

3C.R7 

35.26 

40.00 

25.62 

23.63 

25.67 

2  6.7  r 

26.0  4 

26.60 

27.71 

29.63 

3  2.30 

36.00 

40.00 

25.32 

29.32 

24 ,56 

2  r.44 

24.6  3 

53.04 

30.97 

32.30 

34.28 

37.03 

40.00 

34.46 

34.46 

64.46 

34.52 

34.62 

54.83 

35.26 

36.00 

37.03 

38.46 

40.00 

40.00 

40.00 

40.00 

4  1.00 

43.00 

40.90 

40.00 

40.00 

40.00 

40.00 

40.00 

Sample  of  isotherm  locations  output  from  ADI 
for  rectangle  problem  (file  POINTS) 


3*00 

THE  FOLLOWIMG  144  POINTS  REPRESENT 
90.00  DEGREE  ISOTHERM: 
0.000  2.300 

0.100  2.000 

0.200  2.000 

0 • J  0  0  2.000 

0.400  2.000 


0.500 
0.600 
0.  700 
0.600 
0.900 
1.000 
1.000 
1.000 
1.000 
1.000 
1.000 
1.000 
1.000 
1.000 
1.000 
1.000 
1.000 
1.000 
1.000 
1.030 
1.000 
1.000 
1.000 
1.000 
1.000 


2.000 

2.000 

2.000 

2.000 

2.000 

0.100 

0.200 

0.300 

0.400 

0.500 

0.600 

0.700 

0.800 

0.900 

1.000 

1.100 

1.200 

1.300 

1.400 

1.500 

1.600 

1.700 

1.800 

1.900 

2.000 


TIME  = 


5.00  hours: 


APPENDIX  B:  PROGRAM  ADIPC  AND  SAMPLE  INPUT  AND  OUTPUT 


Program  ADIPC  and  its  subroutines 


1  > 

2)  c  ••***••«••****•****•.«******»•***•*••***•»*»**.•*•..**•**..**..».**••**• 

3 >  C  **»****»«*»•***♦*•»*»•♦**•****«*•******»*******«**•********•*»*«***•»•*» 

4)  C  ADIPC 

5)  C 

6)  C  THIS  FORTRAN  PROGRAM  SOLVES  FOR  TWO  DIMENSIONAL  TEMPERATURE 

7)  C  DISTRIBUTION  RESULTING  FROM  CONDUCTION  HEAT  TRANSFER.  BOUNDARY  CONDITIONS 
H)  C  MAY  INCLUDE  SPECIFIED  TEMPERATURES*  CONVECTIVE  SURFACES*  AND 

9)  C  SEMI-INFINITE  BOUNDARIES.  ADIPC  USES  AN  IMPLICIT  ALTERNATING  DIRECTION 

10)  C  FINITE  DIFFERENCE  NUMERICAL  TECHNIQUE  ON  A  GRID  WITH  SQUARE  ELEMENTS. 

11)  C  THE  MATERIAL  PROPERTIES  <TK*CP.RO*H>  MAY  BE  DIFFERENT  FOR  EACH 

12)  C  NODE  OF  THE  GRID.  AND  MAY  CHANGE  WITH  TIME. 

13)  C  ADIPC  USES  THE  APPARENT  HEAT  CAPACITY  METHOD  FOR  PHASE  CHANGE. 

14)  C  THIS  PROGRAM  WAS  WRITTEN  BY  MARY  REMLEY  ALBERT  AT  CRREL.  1981. 

15) 

16)  C  DATA  FOR  ADIPC  IS  GATHERED  BY  ADDATA*  WHICH  PUTS  IT  INTO  FILE  ADPDAT. 

17)  C  SEE  ADDATA  FOR  AN  EXPLANATION  OF  VARIABLES*  AND  TO  SET  UP  THE  PROGRAM 

18)  C  FOR  YOUR  PROBLEM. 

19)  C  ADJUSTMENTS  FOR  CHANGING  MATERIAL  PROPERTIES  WITH  TIME  (DURING  THE  RUN) 

20)  C  SHOULD  BY  TYPED  INTO  THE  MAIN  PROGRAM  AS  INDICATED  NEAR  LINE  170. 

21)  C 

22)  C  ....***»**»»•»****»**.»****•*******.***********.**.•*•***.**«*»•*.*• 

23)  C 

24)  C 

25)  C 

26)  C  DO  NOT  CHANGE  THE  NEXT  130  LINES. 

27)  IMPLICIT  INTEGER* 2 ( A , 8 ♦ I -L , N *Q , W , X , Y , Z > 

28)  IMPLICIT  DOUBLE  P R EC  I S I  ON (C -H»M • 0  * P « R- V > 

29)  C  Ml  IS  THE  COMMON  LOCATION  FOR  FOR  MAIN  AND  ADDATA  VARIABLES 

30)  C  M12  IS  THE  COMMON  LOCATION  FOR  FOR  MAIN*  ADD A  TA » AND  ISOTHM  VARIABLES 

31)  C0MM0N/M12I/  X*Y*N1S0«Q 

32)  C0MM0N/M12R/  R A Y ( 8 0 *8 0  * 3 ) ♦ T 1  SO (9 > » DS . DEL T 

33)  COMMON/Mll/  IMAX*A,ITRT*KRNR(4)*ITPC 

34)  COMMON/MIR/  DI, TDEL«H( 2) ,FLXT (80 ) . 

35)  &  FLXB (80) ,FLXL (80  )  *FLXR(80) ,TMPT (80)« 

36)  &  TMPB(80).TMPL(80).TMPR(80)*ROPC(2)»CPPC(2)»HL(2>.TPC(2) 

37)  COMMON/MMM /  T EMP ( 80 ,80 > « E A ( 2 0 0 ) * 

38)  &  EB<200),C(200),D(200).CP<«0,80>,TK(80»80), 

39)  &  RO (8  0,80  ) *0LDT (8  0,80 )  ,CPO (8  0,80  )  ,ROO  (80 .80  )  * 

40)  &  IFLAC(10),ISTAT(80,80) 

41)  C  OPEN  FILES 

42)  CALL  C0NTRL(3*'A0PDAT»,5> 

43)  CALL  CONTRL (2,*A0P0UT  * ,6) 

44)  C  FUN1TS  7  «,  8  USED  IN  ISOTHM.  FUNIT  10  NOT  USABLE 

45)  CALL  C0NTRL(2,*ADPNDT» ,9) 

46)  CALL  C0NTRL(2,*ADPTMP»,11) 

47)  CALL  C0NTRL(2,'P0INT1»,13> 

48)  WRITE(l.l) 

49)  1  F0RMAT(1X,»IF  YOU  CHANGEO  ADDATA  SINCE  YOU  LAST  RAN  THIS*** 

50)  C  •  TYPE  "1",»/»1X**  AND  ADDATA  WILL  BE  EXECUTED .» ,/»lX» 

51)  C  'OTHERWISE,  TYPE  "2"*) 

52)  READ(1,»>  IT 

53)  IF(IT.EQ.2>  GO  TO  4 

54)  CALL  ADDATA 

55)  GO  TO  12 

56)  4  CONTINUE 

57)  C  READ  DATA  INPUT  FROM  FILE  ADPDAT 

58)  READ(5,10)  DS .DEL T »DI * T DEL 

59)  10  FORMAT (1X.4F10.5) 

6  G )  READ(5,20)  A ,X , Y , N  ISO , I TRT *  I M AX , I T PC 

61)  20  F0RMAT(1X*7I5> 

62)  R  EAD (5*16)  ( (ROPC (K ) ,CPP C (K ) ,HL ( K > »TP C (K > ) »K =1  * A ) 

63)  16  FORMAT ( 1X.4F10.5) 

64)  READ(5*15)(TIS0(B)*B=1*NIS0) 

65)  15  F0RMAT(1X*F7.2> 

66)  READ(5,30>  ( H(L ) , L  =  1  *  A  ) 

67)  30  FORMAT( 1X.F10.5) 

68)  READ (5*35)  ( ( R A Y ( I  * J* 3 > , J= 1  *  X i , I  =  1  * Y) 

69)  READ ( 5,35 )  ( (RAY ( I*J*2)*J=1*X)*I=1*Y> 

70)  READ(5»35>  ((RAY(I*J*1)*J=1*X)*I=1*Y) 

71)  35  F0RMAT(1X»17F7.2) 

72)  READ(5,34)  (KRNR ( J ) ,J= 1 *4 ) 

73)  34  F0RMAT(1X,4I2) 

74)  READ(5*35)  ( FLX T ( J ) , J= 1 . X > 

75)  READ(5,35>  (FL XB ( J ) * J= 1  *  X ) 

76)  R  EA  D ( 5 , 35  )  ( FL XL ( I ) *  I  =  1 , Y  ) 

77)  READ(5,35>  ( FL XR ( I ) *  I  =  1  * Y > 

78)  READ(5,35)  ( TMP T ( J ) , J= 1 , X ) 

79)  READ(5*35>  ( TMPB ( J ) * J= 1  *  X ) 

80)  READ(5,35)  ( TMPL ( I ) *  I  =  1  * Y ) 


READ  <  5*  35 )  C  TMPR <I),I=1,Y) 

CONTINUE 

WRITE  INITIAL  DATA  INTO  OUTPUT  FILE  ADPOUT 
WRITE  <6,21) 

F0RMAT</»1X, »DATA  FOR  THIS  RUN  OF  AOIPC:*,/,/) 

UR  I TE  <6,22 ) 

FORMAT  </, IX, *  DS  DELT  DI  TDEL*) 

WRITE<6,10>  <OS,OELT,DI,TOEL) 

WRITE  <6,23 ) 

FORMAT!/. 5X.*A*,4X,*X* ,4X , « Y • , 1 X , *N ISO  ITRT  IMAX  ITPC  •» 
WRITE<6,2Q)(A,X.Y,NIS0,ITRT,IMAX«ITPC) 

WRITE  <6,29 ) 

FORMAT!/, IX.*  R  OPC  <K )  CPPC<K)  HL  <K  )  TPC <K > • ) 

WRITE <6,16)  <<ROPC<K) , CPPC < K > »HL < K > ,TPC<K) ),K=1,A) 
WRITE<6,24> 

FORMAT </, IX, *T ISO <B),B=l,NISO:*) 

WR  ITE< 6,15 ) <TISO<B),B=l, NISO) 

WR1TE<6«25) 

FORMAT!/, 3X,»H<K) •) 

WRITE<6,30)  <H<L),L=1,A) 

WRITE  <6, 33)  <KRNR<J),J  =  1,4) 

FORMAT!/, 1X,*KRNR tl)  =  » ,12, »  KRNR < 2 ) =» , 12 , *  KR  NR  <  3)  =  •  » 1 2  , 

&  •  KRNR  <4>=*,I2) 

URI TE  <6 ,26  > 

FORMAT < /,3X, *FLXT  < J) • , 3X , *FLXB < J ) * . 3X » *TMPT < J ) »,3X, 

&  *TMPB<J)*> 

WRI TE <6,2  7) <  <FLXT  < J),FLXB< J) ,TMPT < J ) ,TMPB<J ) ) , J  =  1 ,X ) 
FORMAT<1X,4F10.2) 

WR  I  TE  <6,28  > 

FORMAT!/, 3X,*FLXL<I ) • , 3X , *FL XR < I ) • , 3X , • TMPL < I ) • ,3 X, 

&  • TMPR  < I  )  *  ) 

WRITE <6, 27)  <<FLXL<I)*FLXR<I ) .TMPL < 1 ) ,TMPR ( I ) ) , I  =  1, V  ) 

RRX=X/17 

JJX=RRX 

JJR=X-<17*JJX> 

DO  63  K  =  l,3 

60  TO  <64, 65, 66), K 

WRITE  <6,31) 

GO  TO  67 
URITE<6«45) 

GO  TO  67 
WRITE <6,44) 

DO  66  JJ=1,JJX 
IF(JJX.EQ.D)  GO  TO  6B 
J2  =  l  7  * J J 
J  1  =  J2-16 

-RITE  <6.500)  <<RAY<I, J ,K ) , J= Jl » J2 > , I = 1 , Y ) 

I FIJJ.LT.JJX)  WRITE <6, 503) 

CONTINUE 

J2=X 

Jl=X-UJR+l 

IF(JJR.EQ.O)  GO  TO  69 
IF<JJX.NE.C)  UkITE(6,503) 

DO  69  1=1, Y 

-RITE  <6,500)  <R  AY <  1  ,  J • K )  ,0  =  J 1 , J2 ) 

CONTINUE 

CONTINUE 

FORMAT </,lX, *RA Y< I , J,3 >  NODAL  MATERIAL  TYPE*) 
F0KMAT</,lX,'RAY<I,J,2)  NODAL  LOCATION  TYPE*) 

FORMAT!/, 1X,*RAY< I, J.l)  TEMPERATURES  AT  THE  START!*) 
WRITE <1,11)  ITRT.ITFC 
WRITE <6, 11)  ITRT. ITPC 

FORMATT/,/, IX, 'TEMPERATURES  WILL  BE  PRINTED  EVERY*. 14, 

&  •  TIME  STEPS. *,/,U,»  ISOTHERMS  WILL  BE  LOCATED  EVERY*, 14, 
4  •  TIME  STEPS.*) 

DO  37  K  =  1 , 9 
I FL AG  <  K ) =0 
CONTINUE 
K  ON  T  =  0 
JFLG=1 
Q  =  0 

I TCC=C 
I TPPsO 

D  1 2  =2  .DO  *DI -1 .DO 
K0X  =  1 


NODAL  MATERIAL  TYPE*) 

NODAL  LOCATION  TYPE*) 
TEMPERATURES  AT  THE  START!*) 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
XXXXXXXXXXX  DON*  T  CHANGE  ANYTHING  ABOVE  THIS  LINE.  XXXXXXXXXXXXXXXXXX 


INITIALIZE  SPECIFIC  HEAT  AND  CONDUCTIVITY 
IF  TEMPERATURE  DEPENDENT,  INSERT  EQUATIONS. 
00  46  J=1,X 
DO  47  1=1, Y 
TEMP<I,J)=RAY<I,J,1) 

K  =K  AY  < 1 , J,3  > 

IF<RAY<I,J,1).LT.TPC<K>)  GO  TO  39 
UNFROZEN 
CP<  I,J)=1.160D0 
RO<1,J)=996.20DO 
GO  TO  47 
CONTINUE 

FROZEN 
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174)  CPC  It J)=.5815D0 

175)  R0( It J)=917.00Q 

176)  47  CONTINUE 

177)  46  CONTINUE 

178)  C  •**.»..**•  ****..*.**.  ***********  *********** 

179)  2002  CONTINUE 

180)  C  CONES  HERE  AFTER  EVERY  COMPLETE  TINE  STEP 

181)  C  NOTE  THAT  Q=0  FOR  1ST  ITERATION 

182)  C 

183)  C  IF  SPECIFIED  TEMPERATURES  CHANGE  WITH  TIME*  SPECIFY  THE 

184)  C  CHANGES  HERE.  (NUMBERS  57  AND  58  ARE  FREE  FOR  LOOPS  IF  NEEDED.) 

185)  C 

186)  C 

187)  C  ADJUST  THERMAL  PROPERTIES  WITH  TEMPERATURE. 

188)  C  IF  APPROPRIATE.  USE  RAY(I,J*3)  TO  INDICATE  MATERIAL  TYPE. 

189)  C  (BUT  DO  NOT  CHANGE  THE  NEXT  9  LINES.) 

190)  DO  50  J=1*X 

191)  DO  51  I =1 » Y 

192)  K=RAY(I*J»3> 

193)  CPO(I*J)=CP(I*J) 

194)  ROO(I«J)=ROtI.J) 

195)  CC=TPC(K)*(TDEL/2.D0) 

196)  DD=TPC(K)-(TDEL/2.D0) 

197)  IF(RAY(I*J,l).LT.DO)  GO  TO  60 

198)  IF(RAY( I«J*1) .GT.CC)  GO  TO  61 

199)  C  UNDERGOING  PHASE  CHANGE 

200)  I  ST  AT ( I » J) =3 

201)  TK( I* J)=2.210D0 

202)  GO  TO  51 

203)  60  CONTINUE 

204)  C  FROZEN 

205)  ISTAT(I.J)=1 

206)  TK(I *J)=2.210D0 

207)  RO(I*J)=917.0D0 

208)  CP( I*J)=.58150D0 

209)  GO  TO  51 

210)  bl  CONTINUE 

211)  C  UNFROZEN 

212)  I  ST  AT ( I  * J) =2 

213)  TK( I* J)=.580D0 

214)  RO(I*J)=998.2DO 

215)  CP( I «J>=1 .163000 

216)  51  CONTINUE 

217)  50  CONTINUE 

218)  C 

219)  C 

220)  C  XXXXXXXXXXXXXX  DON • T  CHANGE  ANYTHING  BELOW  THIS  LINE  XXXXXXXXXXXXXXXXXXXX 

221)  C  XXXXXXXXXXXXXX  UNLESS  YOU  KNOW  WHAT  YOU  ARE  DOING.  XXXXXXXXXXXXXXXXXXXX 

222)  2001  CONTINUE 

223)  C  COMES  HERE  AFTER  EACH  PASS 

22 4)  IF(JFLG.EQ.l)  GO  TO  161 

225)  C  FOR  2ND  PASS 

226)  DO  1600  J=1.X 

227)  KONT=C 

228)  DO  41  K=1*Y 

229)  EA(K)=0.0D0 

230)  EB(K)=0.0C0 

231)  C(K)=0.0DG 

232)  D(K)=0.0D0 

233)  D(K)=C.GD0 

234)  41  CONTINUE 

235)  DO  1700  I  =  1*Y 

236)  GO  TO  1601 

237)  161  CONTINUE 

238)  C  FOR  1ST  PASS 

239)  00  160  I  =  1  *  Y 

240)  KONT=0 


241) 

242) 

243) 

244) 

245) 

246) 

247) 

42 

248) 

249) 

1601 

250) 

C 

251) 

252) 

253) 

254) 

255) 

71 

256) 

257) 

258) 

72 

259) 

260) 
261 ) 

73 

262) 

263) 

264) 

74 

265) 

C 

266) 

267) 

C 

DO  42  K=1*X 
EA(K)=0.0D0 
EB(K>=0.000 
C (K ) =G . QDO 
D(K)=0.0D0 
D(K)=O.GDO 
CONTINUE 
DO  170  J=1»X 
CONTINUE 


FIGURE  RESULTANT  CONDUCTIVITIES  BETWEEN  NODES 
KONTrKONT+1 
INDEX=RAY(I«J«3) 

IF(I.EQ.l)  GO  TO  71 

TK1=(2.0D0»TK(1»J)*TK( (I-1)*J))/(TK(1*J)*TK((I-1)»J) ) 
CONTINUE 

IF(J.EQ.l)  GO  TO  72 

TK2=(2.GD0*TK(I»J)*TK( I  »  ( J- 1  )  >  >  /  (  TK  ( I  ,  J  >  *TK  (  I  *  ( J- 1  >  >  > 
CONTINUE 

IF(I.EO.Y)  GO  TO  73 

TK3  =  (2.0D0*TK(I»J)*TK( ( I *1 ) . J ) ) / ( T K ( I , J ) *TK ( ( I ♦ 1 )  ,  J )  ) 
CONTINUE 

IF(J.EQ.X)  GO  TO  74 

TK4  =  (2.000*TK(I*J)*TK(I«(J*1)))/(TK(I*J)4TK(I*(J'»1>)> 
CONTINUE 


I 


FIGURE  R3  FOR  TIME  T-DELT/2*  R4  FOR  TIME  T 
F ( I S T AT (I «J).NE.3)  GO  TO  75 


NOOE  PRESENTLY  PHASE  CHANGE 
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R3  =  2.0C*0S»DS,R0PC(IN0EX>»CPPCUNDEX>/DELT 
R4  =  R3 
GO  TO  78 
75  CONTINUE 

C  NOOE  NOT  PHASE  CHANGE  PRESENTLY 

R3=2.O0*DS*DS«RO0 (I ,J> *CPOU ,J)/DELT 
R4=2.DO»DS»DS*RO< I.J)*CP<I.J)/DELT 
78  CONTINUE 
C  ROUTE  CORNERS 

IF<RAY«I,J,2).GE.5>  GO  TO  38 
I J2=R  A Y ( I , J «  2) 

GO  TO  ( 101,102,103,104) « IJ2 
C  ROUTE  THE  REST  OF  THE  NODES 

38  CONTINUE 

1 J2=RAY<I,J.2)-4 

GO  TO  (105,106,107,108.109,110.111.112,113,114,115, 

A  116.117, 118 ,119, 120.1 21,122 ,123.124, 125 ,126,127, 

&  128.12S, 130, 131,132), IJ2 
C 

110  CONTINUE 
C  SEMI-INFINITE  BOUNDARY 
C  RIGHT  SIDE 

IF(J.NE.X)  GO  TO  1100 
IFLAG(7)=11 

IF(JFLG.NE.l)  GO  TO  1108 
EA(K0NT)=TK2/0I 
EB<K0NT)=-TK2/D1-012«R4 
C(K0NT)=0.D0 

D(KONT)=-TKl*OI2*RAY(( I-1),J»1)-TK3»DI2*RAY( (I ♦ll.J, li¬ 
ft  <TK(I,J)/uI >*TMPR(I),(TK1*DI2*TK3*DI 2-DI 2 *R 3+TK C 1 , J ) /D I ) 
&  *RAY(I«J,1) 

GO  TO  1 b 5 
1108  CONTINUE 

EA(K0NT)=TK1*0I2 

EBIRONT  >=-TKl«D12-TK3»D12-DI2»R4 
C(KONT)=TK3»DI2 

D  (KONT >  =  -<TK2/DI >  *RAY  < I . (J-l > ,1 )-<TK < I ,J >/DI >*TMPR< I )♦ 

&  (TK2/DI*TKU,  J)/01-DI2*R3)*RAYC  I,  J,l> 

GO  TO  165 

1100  CONTINUE 

C  LEFT  SIOE 

IF(J.NE.l)  GO  TO  1101 
I  FLAG (3 ) =1 1 

I F ( JFLG.NE • 1 )  GO  TO  1104 
EA<KONT)=O.DO 
EB(K0NT)=-TK4/0I-DI2*R4 
C(K0NT)=TK4/DI 

DCK0NT)=-TK1*DI2*RAY<< I -1 ) , J ,1 ) -TK 3*D I 2»R A Y ( Cl ♦1),J,1>- 

5  (TK(I,J)/DI)*TMPL(I)*<TK1*DI2*TK3*DI2“DI2»R3»TK( I, J)/DI) 
H  »RAY ( I • J  ,  1  ) 

GO  TO  1 b5 
1104  CONTINUE 

EA(K0NT)=TK1»0I2 

EB(K0NT)=-TK1*DI2-TK3*DI2-D12,R4 

C(K0NT)=TK3*DI2 

0<K0NT)=-(TK4/0I)*RAY< I , < J*1 ) » 1 ) -C TK ( J , J) /DI ) * TMPLI I > ♦ 

S.  <TK4/DI*TK  (I  ,J)/DI-DI2»R3)*RAYU,J*1> 

GO  TO  165 

1101  CONTINUE 

C  TOP  OF  GRID 

IF  (I.NE.l)  GO  TO  1102 

I  FLAG (1 ) =1 1 
URITE(1,9997)  I,J 
GO  TO  9999 

1102  CONTINUE 

C  BOTTOM  OF  GRID 

IFLAG(5)  =  U 

IF< JFLG.NE. 1)  GO  TO  1106 
EA(K0NT)=TK2*0I2 

EB(K0NT)=-TK2*DI2-TK4»DI2-R4*DI2 

CCK0NT)=TK4*DI2 

D(K0NT)--(TK1/DI) *R A Y t < I -1 ) , J, 1 > -< TK ( I , J > /D I ) * TMPBC J ) ♦ 
i  <TK1/DI*TK(I,J>/DI-R3*DI2)*RAY(I,J,1> 

GO  TO  165 
1106  CONTINUE 

EA(K0NT)=TK1/DI 

EB(KONT)=-TK1/OI-R4*DI2 

C(KONT)=O.DO 

D  (K0NT)=-TK2*0I2,RAY<I,«J-1),1)-TK4*0I2*RAY<I,1J4  1),1)- 

6  «TKCI,J)/DI >»TMPB<J)*CTK2*DI24TK4*DI2-DI2*R34TKCI,J>/BT ) 
t  *R AY ( I , J, 1 ) 

GO  TO  165 
109  CONTINUE 

C  CONVECTIVE  SURFACE  BOUNDARY 
C  TOP  CONVECTIVE 

IFLAG!1>=12 

IF< JFLG.NE. 1)  GO  TO  1092 
EA(KONT)=.500*TK2 

EB(KONT>=-.5DO*TK2-.5DO*TK4-.5DO«R4 

C(K0NT)=.5D0»TK4 

0(K0NT)=-TK3*RAY( C I *1 ) , J , 1 ) ♦ I TK 3*H (I NDE X ) «DS-. 5D0 *R3 ) * 

&  RAY<I,J,1>-H<INDEX>*DS*TMPT (J) 

GO  TO  165 
1092  CONTINUE 


-  4  .  *  *-„***’•  S*  v‘ 

'  •  '  •  V  .*  .•  /  .“•*  .*  ■ 
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362) 

E  A(KONT  >  =  0*00 

363) 

EB(KONT)=-TK3-.5DO*R4 

364) 

CtK0NT)=TK3 

365) 

D(KONT)=-.5DO»TK2*RAY(I*(J-1>»1)“.5DO*TK4*RAY{ I,< J*1 >»1>* 

366) 

&  (.5DO*TK2*.5DO*TK4-.5DO*R3*H(INDEX)»DS)*RAY(I *J*1)-H (INDEX) *DS 

36  7 ) 

&  TMPT(J) 

368) 

GO  TO  165 

369) 

132 

CONTINUE 

370) 

C 

RIGHT  SIDE  CONVECTIVE 

371) 

IFLAG(7)=12 

372) 

IF(JFLG.NE.l)  GO  TO  1094 

373) 

EA<K0NT)=TK2 

374) 

EB(K0NT)=-TK2-.5D0*R4 

375) 

C (KONT  >  =0 .DO 

376) 

D(K0NT)=-.5D0*TK1*RAY( (1-1). J. 1 > -.500  * TK3* RA Y< ( 1*1), J,l>* 

377) 

&  (.5DO*TK1*.5DO*TK3-.5DO*R3*H(INDEX)*DS)*RAY(I,J»1>- 

378) 

&  H(INDEX)*DS*THPR(I) 

379) 

GO  TO  165 

38  0) 

1094 

CONTINUE 

381 ) 

EA(KONT)=.5DO*TK1 

382) 

EB<KONT)=-.5D0*TKl-.5O0»TK3*.5D0*R4 

383) 

C  (KONT)  =  .5DO*TK3 

384) 

D(KONT)=-TK2*RAY(I»(J-l)»l)*tTK2-.5D0*R3*H (INDEX) *DS>* 

385) 

&  RAY( I»J.1)-H(INDEX)*DS*TMPR (I) 

366) 

GO  TO  165 

387) 

131 

CONTINUE 

368) 

C 

BOTTOM  CONVECTIVE 

389) 

IFLAG(5)=12 

390 

IF( JFLG.NE.l)  GO  TO  1096 

391) 

EA(KONT)=.5DO*TK2 

392) 

EB(KONT)=-.5D0*TK2-.5O0*TK4-.5D0*R4 

393) 

C(KONT)=.5D0*TK4 

394) 

D(K0NT)=-TK1*RAY((I-1) • J. 1 )♦( TK 1 *H ( INDEX ) *DS-. 5D0 *R3  )  * 

395) 

t  RA Y( I  * J.l ) -H(INDEX)*DS*TMPB (J ) 

396) 

GO  TO  165 

397) 

1096 

CONTINUE 

398) 

EA(K0NT)=TK1 

399) 

EB(KONT)=-TK1-.5DO*R4 

400) 

C ( KONT ) =0 • 00 

401) 

D(KONT)=-.5DO*TK2*RAY(It(J-l)*l)-.5DO*TK4*RAY(I.(J*l)»l)4 

402) 

g  (•5D0*TK2*.5D0*TK4-.5D0*R3*H(INDEX)*DS)*RAY(I»J»1)-H(INDEX)*DS 

403) 

fc  TMPB ( J  > 

404) 

GO  TO  165 

405) 

130 

CONTINUE 

406) 

C 

LEFT  SIDE  CONVECTIVE 

407) 

IFLAG(3)=12 

408) 

IF( JFLG.NE.l)  GO  TO  1097 

409) 

EA(KONT)=O.DO 

410) 

EB(KONT)=-TK4-.500*R4 

411) 

C(K0NT)=TK4 

412) 

D(KONT)=-.5O0*TKl*RAY( ( I -1 ) « J » 1 ) -.5D0 *TK 3*R AY ( (1*1). J.l)* 

413) 

&  (.5DO*TK1*.5DO*TK3-.5DO»R3*H(INDEX)*DS)*RAY(I, J.l)- 

414) 

&  H< INDEX) *DS *T HPL ( I ) 

415) 

GO  TO  165 

416) 

1097 

CONTINUE 

417) 

EA(KONT)=.5DO*TK1 

418) 

EB(KONT)=-.500*TK1-.500*TK3*. 5D0*R4 

419) 

C(KONT)=.500*TK3 

420) 

0  (KONT)--TK4*RA Y( It (J*l) *  1 )♦( TK4-. 5D0*R3*H( INDEX) *DS)* 

421) 

&  RAY( I.J,1)-H(INDEX)*DS*TMPL(I) 

422) 

GO  TO  165 

423) 

107 

CONTINUE 

424) 

C  CONSTANT  TEMP  BOUNDARY 

425) 

EB( KONT ) = 1 

426) 

D(KONT)=RAY(I*J«l) 

427) 

IF(I.NE.l)  GO  TO  1070 

428) 

IFLAG(l)slO 

429) 

GO  TO  165 

430) 

1070 

CONTINUE 

431) 

IF(J.NE.l)  GO  TO  1071 

432) 

IFLAG<3)=10 

433) 

GO  TO  165 

434) 

1071 

CONTINUE 

435) 

IF(I.NE.Y)  GO  TO  1072 

436) 

IFLAG(5)=10 

437) 

GO  TO  165 

438) 

1072 

CONTINUE 

439) 

IF(J.NE.X)  GO  TO  1073 

440) 

IFLAG(7)=10 

441) 

GO  TO  165 

442) 

1073 

CONTINUE 

443) 

106 

CONTINUE 

444) 

C  CONSTANT  INTERIOR  NODE 

445) 

EB(KONT)=l 

446) 

D(KONT)=RAY(I«J*l) 

447) 

GO  TO  165 

448) 

105 

CONTINUE 

449) 

C  VARIABLE  INTERIOR  NODE 

450) 

I F ( JFL6.NE .1 )  GO  TO  1050 

451) 

EA(K0NT)=TK2 

452) 

C (KONT ) =TK4 

453) 

EB(K0NT)s-TK2-TK4-R4 

454) 

D(K0NT)=-TK1*RAY( (1-1) » J.l ) -TK3*RA Y ( ( I *1 1 , J* l)  ♦ 

455) 

t  (TK1*TK3-R3)*RAY(I*J»1) 

456)  GO  TO  165 

457)  1 05  G  CONTINUE 

458)  EA(K0NT)=TK1 

459)  C(K0NT)=TK3 

460)  EB(K0NT)=-TK1-TK3-R4 

461)  0<K0NT)=-TK2»RAY<  I  « < J- 1 ) « 1 ) -T K 4* R A Y ( 1 ♦ ( J* 1 ) *  1 ) ♦ < T K2*TK4-R 3 ) * 

462)  &  R  A  Y  (  I,  J,l) 

463)  GO  TO  165 

464)  129  CONTINUE 

465)  C  CONSTANT  HEAT  FLUX 

466)  C  RIGHT  SIDE 

467)  IFLAG<7)=13 

468)  I F< JFLG.NE . 1 )  GO  TO  1301 

469)  E  A  (  R  0  N  T )=TK2 

470)  EBfKONT > = -T K 2-. 50 0 *R4 

471)  0  (KONT)=-.5DO*TK1  *RAY(  ( I -1 > ♦ J* 1 > -.500 *TK 3*R AY(  (1*1) *0*1)* 

472)  &  <.5DO*TK1*.5DO»TK3-.5DO*R3)*RAY(I»J*1)-FLXR(I)*DS 

473)  C(KONT)=O.DO 

474)  GO  TO  165 

475)  1301  CONTINUE 

476)  EA(K0NT)=.5DQ«TK1 

477)  EB(KONT)=-.5DO*TK1-.5DO*TK3-.5DO»R4 

478)  C (KONT)=.5D0«TK3 

479)  D(K0NT)=-TK2*RAY(I*(J-1)»1)*CTK2-.5D0*R3)*RAY(I*J*1) 

480)  &  -FLXK  < I ) *  DS 


481) 

482) 

483)  C 

484) 

485) 

486) 

487) 
488  ) 
4b9) 

49  0) 

491) 

492) 


GO  TO  1»5 

127  CONTINUE 

LEFT  SIDE  CONST  FLUX 
IFLAG<3)=13 

IF(JFLG.NE.l)  GO  TO  1303 
EA(KONT>=0.00 
Eb(K0NT)=-TK4-.50G*R4 
C(K0NT)=TK4 

D(KONT)=-.bDOMKl»KAY(  (l-l).J,l)-.bD0*TK3*RAYC  (1*1),  J,l)* 

6  <.5D0*TKl*.5O0*TK3-.5D0*R3)*RAY(I*J»l>-FLXL(I>*DS 
GO  TO  165 
1303  CONTINUE 

EA(KONT)=.bDO*TKl 

EB(KONT)=-.5DO*TK1-.5DO*TK3-.5DO*R4 

C(KONT)=.5DO*TK3 

0<K0NT)=-TK4*RAY( I, ( J* 1 ) ,1 )  ♦  C  T K4 -  .  5 DO  *R 3  >  * R  A  Y  ( I . J , 1 ) -FL XL C I  )»DS 
GO  TO  165 
108  CONTINUE 

TOP  CONST  FLUX 
IFLAGU  )  =  13 

IF(UFLG.NE.l)  GO  TO  1305 
EAIKONT )=.5D0*TK2 

EB(KONT)=-.5D0*TK2-.5D0*TK4-.5D0*R4 

C(K0NT)=.50C*TK4 

D(KONT)=-TK3*RAY( (1*1) ♦ J * t ) ♦ ( TK 3-. 5D0 *R 3 > *R A Y ( I » J ,1 > -FL XT C J ) *DS 
GO  TO  165 
1305  CONTINUE 

EA(nONT)=O.DQ 
EBCKONT )=-TK3-.5D0*R4 
CIK0NT)=TK3 

0(K0NT)=-.5D0*TK2*RAY(I«(J-l)*l)-.5D0*TK4*RAY(I*(O*l)«l)* 
l  < .500*TK2*.500*TK4-.50G*R3)*RAY  (I* J«1)-FLXT CJ  )*DS 
GO  TO  165 

128  CONTINUE 

BOTTOM  CONST  FLUX 
I F( JFLG.NE. 1 )  GO  TO  1307 
IFLAG(5)=13 
E  A  <  K ONT )=.5DC*TK2 

EBIKONT >=-.5D0*TK2-.5O0*TK4-.5D0*R4 
C(KONT)=.5DO*TK4 

0(KONT)=-TKl*RAY( II -1 ) * J . 1 ) ♦ ( TK 1 -. 5D0*R 3 ) *R A Y( 1*0*1) 

A  -FLXB(J)»DS 
GO  TO  165 
1307  CONTINUE 

EA(K0NT)=TK1 

EB(KONT)=-TK1-.5DO*R4 

C(KONT)=O.DO 

0<KONT)--.5DO»TK2*RAY( I  *  I J-l > , 1 ) -.500 *TK 4*R AY (I  *  1 0*1 ) *1 ) ♦ 

&  (.5D0»TK2*.5D0»TK4-.5D0*R3)*RAY(I*J«1> 

4  -FLX8 ( J  >  *DS 
GO  TO  165 

111  CONTINUE 

NODE  ADJ  TO  LEFT  SEMI-INF 
I F (JFLG.NE • 1 )  GO  TO  1110 
EA(K0NT)=TK2/0I 
EBUONT  >=-TK4-TK2/DI-R4 
C(K0NT)=TK4 

D(K0NT)=-TK1*RAY((I-1)*J«1)-TK3*RAY<(I*1)»J*1)*(TK1* 
l  TK3-R3)»RAY(I»J»1) 

GO  TO  165 
1110  CONTINUE 

EA(K0NT)=TK1 

EB(K0NT>=-TK1-TK3-R4 

C(K0NT)=TK3 

D(KONT)=-TK4*RAYC I  * < 0* 1 > « 1 > - ( TK2/0 1 > *R AY < I  * ( 0- 1 > *  1 ) * 
t  (TK4*TK2/DI-R3)»RAY(I « J « 1 ) 

GO  TO  165 

112  CONTINUE 

NODE  ADJ  TO  BOTTOM  SEMI-INF 
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550)  1 F  < JFLG  .NE  •  1 )  GO  TO  1120 

551)  £  A  < KONT >  =  TK2 

552)  EBtKONT > =-TK 2-T K4 -R 4 

553)  C  t  K  ONT )  =  TK  4 

55 A)  D tKONT)=-TKl*RAYt  1 1-1 >»J.l)-tTK3/DI >*RAYttI*l>, J.l»* 

555)  &  <TK1*TK3/DI-R3)*RAY(I *  J  .  1 ) 

556)  GO  TO  165 

557)  1120  CONTINUE 

558)  £A(K0NT)=TK1 

559)  EBtKONT >=-TKl-TK3/0I-R 4 

560)  C(KONT)=TK3/Dl 

561)  DtKONT)=-TK2*RAY II  *  tj-l) »1 )-TK4*RAY  t I . t J*1 > ♦ 1 » *f TK2* 

562)  &  TK4-R3)*RAYt I , J,  1) 

563)  GO  TO  165 

564)  113  CONTINUE 

565)  C  NODE  ADU.  TO  RIGHT  SEMI-INF 

566)  IFt JFLG.NE. 1 >  GO  TO  113C 

567)  E  A  t  KONT )  =  TK  2 

568)  EBtKONT )=-TK2-TK4/GI-R4 

5b9 )  C(K0NT)=TK4/DI 

57  0)  D tKONT)=-TKl*RAYt  II -1) , J.  1 ) -TK3*R A Y 1 1 1*1 > » J.l)* 

571)  &  (TK1+TK3-R3)*RAY ( 1  *  J  *  1 ) 

572)  GO  TO  165 

573)  1130  CONTINUE 

574)  EAIKONT )  =  TK  1 

575)  EBtKONT )=-TKl-TK3-R4 

576)  CIKONT)=TK? 

577)  0tK0NT)=t-TK4/DI)*RAYt 1 . t J* 1 > . 1 ) -T K2*R AY t I . t J- 1 > , 1 > ♦ 

576)  &  ITK2*TK4/DI-R3)*RAYC I  *  J » 1 ) 

579)  GO  TO  165 

580)  115  CONTINUE 

581)  C  SQUARE  NODE  ADJ  TO  2  SEMI-INF  SIDES 

562)  IFtJ.NE.2)  GO  TO  1150 

583)  IFtI.NE.2)  GO  TO  1151 

584)  C  TOP  LEFT 

585)  UR  I TE  < 1  *  9997 )  I.J 

566)  GO  TO  9995 

587)  1151  CONTINUE 

588)  C  BOTTOM  LEFT 

589)  IF(JFLG.NE.l)  GO  TO  1152 

590)  EACK0NT)=TK2/DI 

591)  EBIKONT )=-TK2/DI-TK4-R4 

592)  C<KGNT)=TK4 

593)  0<K0NT)=-TKl.RAYt<I-l>,J,l)*tTK3/Dl)*RAYttI*l),J,l)*tTKl* 

594)  &  TK3/DI-R3)*RAYC I.J.l) 

595)  GO  TO  lb5 

596)  1152  CONTINUE 

597)  EACK0NT)=TK1 

598)  EB(KONT)=-TKl-TK3/DI-R4 

599)  C<KONT)=TK3/DI 

60  0)  DCK0NT)=-(TK2/DI ) *R AY t I » C J-l ) « 1 ) -TK4*R AY ( I « C J* 1 )  . 1 ) ♦ 

601)  &  tTK2/01*TK4-R3)*RAYtI. J.l) 

602)  GO  TO  165 

603)  1150  CONTINUE 

604)  IFCI.NE.2)  GO  TO  1153 

605)  C  TOP  RIGHT 

606)  UR  I TE  (  1 .999  7  )  I.J 

607)  GO  TO  9599 

608)  1153  CONTINUE 

609)  C  BOTTOM  RIGHT 

610)  1FCJFLG.NE.1)  GO  TO  1154 

611)  EA(KONT)=TK2 

612)  EBtK0NT)=-TK2-TK4/0I-R4 

613)  C(K0NT)=TK4/0I 

614)  DtKONT)=-TKl*RAYt ( I-l).J.l)«tTK3/DI)*RAYt(l*l)«J. 1)* 

615)  &  tTKl*TK3/DI-R3)*RAYtI » J . 1 ) 

616)  GO  TO  165 

617)  1154  CONTINUE 

618)  EAtKONT ) =TK1 

619)  EBtK0NT)=-TKl-TK3/DI-R4 

620)  CtKONT)=TK3/DI 

621)  D(KONT)=-TK2*RAY(  I.  (J-l)  ,1)-(TK4/DI  )*RAY(I.(  J*1>,1>* 

622)  &  tTK2*TK4/DI-R3)*RAYtI,J,l) 

623)  GO  TO  165 

624)  121  CONTINUE 

625)  C  RIGHT  SIDE  CONST  FLUX  ADJ.  TO  BOTOM  SEMI-INF 

626)  IFtJFLG.NE.l)  GO  TO  1210 

627)  EAtK0NT)=TK2 

628)  EBtKONT )=-TK2-.5D0*R4 

629)  CtKONT)=O.DO 

630)  D(K0NT)=-.5D0*TKl*RAYt  t I -1 ) . J. 1 ) -< TK3/<2 .00*01) ) *R A Y t t I *1 ) « J. 1 ) - 

631)  A  FLXR  tl )*DS*t.5D0*TKl*TK3/t2.D0*DI )-.500*R3)*RAYt 1 . J . 1 ) 

632)  GO  TO  165 

633)  1210  CONTINUE 

634)  EA(KONT)=.5DO*TK1 

635)  EBtKONT >=-.5O0*TK 1-TK 3 / t 2 .00 *DI ) -.500 *R4 

636)  CtKONT)=TK3/t2.DO*DI> 

637)  DtK0NT)=-TK2*RAYt I « t J- 1 ) « 1 ) -FLXR t I ) *0S* t TK2-.5D0*R3 ) *R A Yt I , J. 1 ) 

638)  GO  TO  165 

639)  123  CONTINUE 

640)  C  LEFT  SIOE  CONST  FLUX  ADJ  TO  80‘TTOM  SEMI-INF 
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641)  IF(UFLG.NE.l)  GO  TO  1230 

642)  EA<KONT)=O.DO 

643)  EB(KONT)=-TK4-.5D0«R4 

644)  C<K0NT)=TK4 

645)  D<KONT)=-.5DO*TK1*RAV< < I - 1 ) , J, 1 ) - ( T K3 / < 2 . D 0  * DI ) ) «RAY ((I*1)«J«1) 

646)  &  FLXL<I)*DS*C.500*TK1*TK3/I2. 00*D1)-. 500 *R3> *R AY( It J» 1 > 

647)  GO  TO  165 

648)  1230  CONTINUE 

649)  EAIKONT)=.5DO*TK1 

650)  EB(KONT)=-.5DO*TK1-TK3/<2.DO*DI)-.5DO*R4 

651)  C(KONT)=TK3/(2.DO*OI> 

652)  D(K0NT)=-TK4*RAY< It <J*1) tl)-FLXL<I )*DS*(TK4-.5D0*R3)»RAY(I«J*1) 

653)  GO  TO  165 

654)  116  CONTINUE 

655)  C  SEMI-INF  NODE  ABOVE  SEMI-INF  CORNER  NODE 

656)  IF(J.NE.l)  GO  TO  1165 

657)  C  BOTTOM  LEFT 

658)  IF(JFLG.NE.l)  GO  TO  1160 

659)  EA(KONT)=O.ODO 

660)  EB(K0NT)=-TK(ItJ)/DI-TK4/DI*R4*DI2 

661)  C(K0NT)=TK4/DI 

662)  D<K0NT)=-TK1*DI2*RAY< ( I -1 ) . J ♦ 1 > - ( TK 3*012 /D I )*RAY( <I*l)t J»l) 

663)  &  -ITK (I tJ)/DI)*TMPL<I )  ♦ <  TK  1  * DI  2*TK 3 *DI2/DI -R 3*DI 2 ) *R AY C 1 1 J 1 1  > 

664)  GO  TO  165 

665)  1160  CONTINUE 

666)  EA<K0NT)=TK1*DI2 

667)  EBCK0NT)=-<TKl*TK3/DI-R4)*DI2 

668)  C(K0NT)=TK3*DI2/DI 

669)  D(K0NT)=-(TK(ItJ)/DI)*TMPLCI)-(TK4/DI)*RAY<Iti J*l)tl)» 

670)  &  <TK(I«J)/DI+TK4/0I-R3*DI2)*RAY<ItJ»l> 

671)  GO  TO  165 

672)  1165  CONTINUE 

673)  C  BOTTOM  RIGHT 

674)  I F(JFLG.NE.l)  GO  TO  1166 

675)  EA(KONT)=TK2/DI 

676)  EBIK0NT)=-TKIItJ)/Dl-TK2/DI*R4*DI2 

677)  C<KONT)=O.ODO 

678)  D<K0NT)=-TKl*DI2*RAY((I-l),J«l)-(TK3«D12/0I)*RAYt(I*l)»J»l) 

679)  &  -<TK<It J)/DI)*TMPL<I)*<TKl*DI2*TK3*DI2/DI-R3*DI2)*RAYCIt Jt  1) 

680)  GO  TO  165 

681)  1166  CONTINUE 

682)  EACK0NT)=TK1*D12 

683)  EB(K0NT)=-<TK1*TK3/DI-R4)*DI2 

684)  CCK0NT)=TK3*DI2/0I 

685)  D(K0NT)=-«TK<ItJ)/0I)*TMPR(I)-CTK2/DI)*RAYCI ti J-l>tl)t 

686)  &  (TK(I»J)/DI*TK2/DI-R3*DI2)*RAY(I*J«1) 

687)  GO  TO  165 

688)  117  CONTINUE 

689)  C  SEMI - INF  NODE  TO  LEFT  OF  SEMI-INF  CORNER  NODE 

690)  IFC JFLG.NE.l)  GO  TO  1170 

691)  EA<K0NT)=TK2*012 

692)  EB(K0NT)=I-TK4/0I-TK2-R4)*DI2 

693)  C(K0NT)=TK4*0I2/D1 

694)  D»K0NT)=<-TK1/DI)*RAYC < I -1 ) t Jt 1 ) -C TK C I t J )/DI )* TMPBC J) 

695)  l  ♦CTKl/DI+TK(I»J)/DI-R3»DI2J*RAY<I»J«i) 

696)  GO  TO  165 

697)  1170  CONTINUE 

698)  EA«KONT)=TKl/OI 

699)  EB(K0NT)=-(TK1*TKCI»J) )/0I-R4*DI2 

700)  C(KONT)=O.ODO 

701)  DtKONT)  =  C-TK4*OI2/DI)*RAYUt  CJ*1  )»  l )-TK2* D 12 »R AY C  I,C  J-D.l)- 

702)  l  (TK<ItJ)/OI)«TNPB<J)*»TK2/OI»TK4-R3)*OI2*RAY(ItUtl) 

703)  GO  TO  165 

704)  119  CONTINUE 

705)  C  SEMI-INF  NODE  TO  RIGHT  OF  SEMI-INF  CORNER  NODE 

706)  I F < JF LG »NE • 1 >  GO  TO  1190 

707)  EACKONT)=TK2*DI2/OI 

708)  EB<K0NT)=<-TK2/0I-TK4-R4)«DI2 

709)  C<K0NT)=TK4*DI2 

710)  DCKONT)  =  t-TKl/DI )*RAY<  < I -1 ) t Jt 1 ) -( TK C I tJ)/DI )*TNPBCJ)* 

711)  &  CTK1/DI*TKCI,J)/DI-R3*DI2)*RAYCI,J»1) 

712)  GO  TO  165 

713)  1190  CONTINUE 

714)  EA<KONT)=TKl/OI 

715)  E8(K0NT)=-(TKl*TK(ItJ) )/DI-R4*DI2 

716)  C<KONT)=O.ODO 

717)  D<K0NT)=C-TK2»DI2/DI)*RAYCI t CJ-1 )  1 1 >-TK4 *DI 2 *R A Y < 1 1 « J*1 ) 1 1 ) - 

718)  t  (TK(I*J)/DI)*TMPB(J)*ITK2/DI*TK4-R3)»DI2*RAY<I«J»1) 

719)  GO  TO  165 

720)  120  CONTINUE 

721)  C  TOP  CONST  FLUX  NODE  ADJ  TO  RIGHT  SEMI-INF 

722)  IF(JFLG.NE.l)  GO  TO  1200 

723)  EA(KONT)=.50DO*TK2 

724)  EB<KONT)=-.50D0*(TK2*TK4/DI*R4) 

725)  C(KONT)=TK4/(2.0DO*DI> 

726)  D<K0NT)=-TK3*RAY< U»l) t J t 1 > -DS*FLXT C J )♦* TK3-.50DO *R3 ) 

727)  t  *R  AY  C 1 1 Jt 1 ) 

728)  GO  TO  165 

729)  1200  CONTINUE 

730)  EA<KONT>=0.000 

731)  EBIKONT)=-TK3-.50D0«R4 

732)  C(K0NT)=TK3 

733)  D(KONT)=-l.*(TK4/(2.0DO*OI))*RAY(l«<J*l>tl>-»5*TK2 
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734) 

735) 

736) 

737) 

738) 

739) 
740  > 

741) 

742) 

743) 

744) 

745) 

746) 

747) 

748) 

749) 

750) 

751) 

752) 

753) 

754) 

755) 

756) 

757) 

758) 

759) 

760) 

761) 

762) 

763) 

764) 

765) 

766) 

767) 

768) 

769) 

770) 

771) 

772) 

773) 

774) 

775) 

776) 

777) 

778) 

779) 

780) 

781) 

782) 

783) 

784) 

785) 

786) 

787) 

788) 

789) 

790) 

791) 

792) 

793) 

794) 

795) 

796) 

797) 

798) 

799) 
800  ) 
801) 
802) 

803) 

804) 

805) 

806) 

807) 

808) 

809) 

810) 
ail) 
812) 

813) 

814) 

815) 

816) 

817) 

818) 

819) 

820) 
821) 
822) 

823) 

824) 

825) 

826) 
827) 


&  *R AY(1 » (J-l ) »1 >♦ (TK4/ (2.0D0*DI >*.50D0*TK2-.50D0*R3> 

5  *RAY(I, U«1)-DS«FLXT( J) 

122  CONTINUE 

C  TOP  CONST  FLUX  AO J  AC  TO  LEFT  SEMI-INF 

IF(JFLG.NE.l)  GO  TO  1220 
EA(KONT)=TK2/(2.CDO*DI > 

E6(KONT>=-.5000»(TK2/DI»TK4*R4) 

C (KONT  >=.50D0*TK4 

0(KONT)=-TK3*RAY ( C I ♦ 1 ) »J«1)-DS*FLXT  <  J )  * 

6  (TK3-.50DO»R3)*RAY(I«J»1> 

GO  TO  165 

1220  CONTINUE 

EA!KONT)=O.ODO 
EB(K0NT)=-TK3-.5»H4 
C(KONT  )=TK3 

0(K0NT)=-(TK2/(2.0D0«D1)>*RAY<1«(J-1)»1)-»50D0*TK4* 

&  RAY< I. Cd*l ).!)♦< TK2/<2. 000*01)*. 5000 *TK4-. 5000*R3) 

&  »RAY(I«J»1)-DS»FLXT(J) 

GO  TO  165 
114  CONTINUE 
118  CONTINUE 

124  CONTINUE 

125  CONTINUE 

126  CONTINUE 

WRITE! 1*9997)  I.J 
GO  TO  9999 

C  ((((((<((((((((((((((((((  ( 

C  EQUATIONS  FOR  THE  CORNERS  OF  THE  GRID 

101  CONTINUE 

C  TOP  LEFT-HAND  CORNER 

GO  TO  (201«202. 203.204,205. 206. 207, 208t209»210)»KRNR«l> 

202  CONTINUE 

C  CONSTANT  FLUX  ON  BOTH  SIDES 

I  FLAG  (2  >  =  1 

IF(JFLG.NE.l)  GO  TO  2020 
E  A ( K  ONT )  =  0  •  D  0 

EB(KONT)=-.5O0«TK4-.25DO*R4 
C (KONT )=.5D0«TK4 

D(KONT)=-.5DO*TK3*RAY( ( I ♦! > , U, 1 ) ♦ ( . 5D C*TK3- ,25D0*R3 ) *R A Y ( I , J , 1 > - 
&  .5D3*0S*(FLXT(J)-»FLXL(I)> 

GO  TO  165 
2020  CONTINUE 

EA(KONT)=O.DO 

EB(KONT)=-.5DO*TK3-.25DG*R4 

C(KONT)=.5DO*TK3 

DCKONT)=-.5D0*TK4*RAY( I , < J*1 ) ♦ 1 ) ♦ t . 500 »TK4-. 2500 *R3) * 

&  RAY<  I,  J,ll-.5O0*OS*<FLXTtJMFLXL(I  )  ) 

GO  TO  165 

203  CONTINUE 

C  CONVECTIVE  ON  BOTH  SIDES 

I FL  AG  <  2 )  =  2 

IFC JFLG.NE.l)  GO  TO  2030 
EA(KONT)=0.00 

EB(KONT)=-.5DO«TK4-.25DO*R4 

C(KONT)=.5DO*TK4 

D(KONT)=-.5DO*TK3*RAY< ( I ♦ 1 > » J. 1 ) ♦ f . 5D 0* TK3*H < I NDE X) *DS- .25D0*R3 ) * 
&  RAY<I.J,1)-.500*H<INDEX)*TMPT<J)-.5DO*H(INDEX>*TMPL<I> 

GO  TO  165 
2030  CONTINUE 

E A (KONT )  =  0  .00 

EB(KONT)=-.5DO*TK3-.25DD*R4 

C(KONT)=.5DO*TK3 

0(KONT)=-.500*TK4*RAY( I  ,  <  J*  1 )  » 1  >  ♦  l  .  5D0»T  K4-.HC  I  NDEX)  *  DS-.25D0  *R3  )  * 
&  RAY( I«J»l)-.5O0*H (INDEX  )*(THPT (J)*TMPLCI)  ) 

GO  TO  165 

204  CONTINUE 

C  SEMI-INF  ON  BOTH  SIDES 

I FL AG ( 2 ) =3 
WRITE <1«9997)  I.J 
GO  TO  9999 
206  CONTINUE 

C  VERT  CONST  FLUX—  HORIZ  CONVECT 

IFL AG(2)=4 

I F ( JFLG.NE.l)  GO  TO  2060 
EA(K0NT)=0.0D0 

EB(KONT)=-.50D0*TK4-.250D0*R4 

C(KONT)=.50D0*TK4 

D(KONT)=-.50D0*TK3*RAY ( ( I ♦ 1 ) , J, 1 ) - . 50D0*FL XL ( I )»DS 
l  -.50 DO »H< INDEX ) * DS *TMPT ( J ) ♦ ( .5000* TK3*. 50 DO^H ( INDEX )*DS 
t  -.2 5000»R3)»RAY( I. J,1 ) 

GO  TO  165 
2060  CONTINUE 

EA(KONT>=0.000 

EB(KONT)=-.5000*TK3-.5000«HCINDEX)*DS 

C(KONT)=.50D0*TK3 

0(KONT)=-.5000*TK4*RAY (I ,  ( J*l).l ) - .50 D0*FL XL ( I )*DS 
&  -.50D0»H( INDEX >»OS*TMPT<J)»(.5000*TK4-.25000*R3> 
t  •RAY(IfJtl) 

GO  TO  165 

205  CONTINUE 

C  VERT  CONVECT— HORIZ  CONST  FLUX 

IFLAG<2>=7 

I F ( JFLG.NE . 1 >  GO  TO  2050 
£A(KONT)=0.000 
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E B(KONT >=-. 5000 *TK4-. 5 000*H ( INDEX >*DS-. 25000 *R 4 
C(KONT)=.50D0«TK4 

D(KONT)=-.50D0*FLXT(J>*DS-.50D0*TK3*RAY( (1*1). J.l)* 
l  (.5000 *TK3-. 2500 0*R3>*RAY(I*J.1)-.50D0*H( INDEX) *DS*THPL( I) 

60  TO  165 
2050  CONTINUE 

E  A<  KONT  >=0.000 

EB(KONT)=-.5QD0*TK3-.25ODO*R4 

C(KONT)=.50D0*TK3 

D(KONT)=-.50D0*TK4»RAY(lt(J*l)tl>-.50D0*H(INDEX>*DS*THPL(I> 

A  *(.50D0*TK4*.50D0*H(INOEX)*DS-«50D0*FLXT(J)*OS-«250O0*R3)* 

A  RAY (  I  « J.l ) 

GO  TO  165 

206  CONTINUE 

:  VERT=SEHI-INF.HORIZ=CONST  FLUX 

IF( JFLG.NE.l)  GO  TO  2080 
EA(K0NT)=0.0D0 

EB(KONT)  =  -(TK ( I «J  >/(2.OD0*DI >*TK 4/ (2 . OOO *01 > 

A  ♦  • 5000*01 2*R4  > 

C(K0NT)=TK4/<2. 000*01) 

D(KONT)=-FLXT(J>*DI2*OS-TK3*DI2*RAY( (I*1>«J«1>* 

(  (TK3*DI2-.5000*OI2*R3>*RAY(I»J»1>-(TK(I .  J  >  /  (2  .  OD  0*0 1 >>*THPLC1 ) 

GO  TO  165 
2060  CONTINUE 

EA(KONT)-O.ODO 

EB(KONT)=-(TK3*.50D0*R4)*DI2 

C(K0NT)=TK3*0I2 

0(K0NT)=-FLXT( J > *01 2«D S- ( TK ( I ♦ J) / ( 2 .0D0*DI > > *T HPL ( I > 

A  -(TK 4/ (2. OOO *01 ) ) *R AY ( I i ( J*1 ) « 1 )♦ ( TK ( I *J ) / ( 2. ODO *DI ) ♦ 

A  TK4/<2. 000*01)-. 50D0*DI2*R3)*RAY< I* J.l) 

GO  TO  165 
210  CONTINUE 

:  VERT=SEHI-INF,HORIZ=CONVECT 

IFLAG<2)=8 

IF! JFLG.NE.l)  GO  TO  2100 
EA(KONT>=0.000 

EB(K0NT)=-(TK4*TK(I.J> )/(2.0D0*DI)-.50D0*DI2*R4 
C ( KONT )sTK4/ (2. 000*01) 

D(KONT>=-H(INDEX)*DI2*THPT(J>-(TK(I»J)/(2.0DO*OI> )*THPL(I) 

A  -TK3*DI2*RAY((I*1)»J»1)*(TK3*DI2*H(INDEX)*0I2-.50D0*DI2*R3) 

A  *RAY(I.J.1> 

GO  TO  165 
2100  CONTINUE 

E  A (KONT )=0 . ODO 

EB(  KONT >=-TK3*DI2-H( INDEX) *DI2-.50D0*DI2*R4 
C(K0NT)=TK3*DI2 

0  (KONT  )  =  -H( INDEX) *D1 2* TMPT ( J >- ( TK4 / (2 .0D0*DI ))*RAY(I ,(J*1> ,1) 

A  -(TK(I,J)/(2.000*D1>>*TNPL(I>*((TK4*TK(I.J>>/(2.0DO*DI>- 
A  .50D0«DI2*R3)*RAY(ItJ»l) 

GO  TO  165 

207  CONTINUE 
209  CONTINUE 

WRITE! 1*9997)  I  *  J 
GO  TO  9999 
201  CONTINUE 

:  CONSTANT  CORNER  NOOE 

EB ( KONT ) =1 
0(K0NT)=RAY( I* J.l ) 

I  FLAG (2 )  =  1 0 
GO  TO  165 
102  CONTINUE 

:  BOTTOM  LEFT-HAND  CORNER 

GO  TO  (211. 212. 213.214. 215.216. 217. 218»219.220>.KRNR(2> 

212  CONTINUE 

:  CONST  FLUX  ON  BOTH  S10ES 

I  FLAG ( 4)  =  1 

IF( JFLG.NE.l)  GO  TO  2120 
E A (KONT ) :Q • 00 

EB(KONT)=-.5DO*TK4-.25DQ*R4 

C(KONT)S,5DO*TK4 

D(KONT)=-.500*TK1*RAY( ( I -1 > » J.l > ♦ ( . 5D0*TK1 -.2500 *R3 > *RA Y ( I . J . 1 > 

A  -.5*0S*(FLX8(J)*FLXL(I>) 

GO  TO  165 
2120  CONTINUE 

EA(K0NT)=.5D0*TK1 

EB(KONT)=-.50Q*TK1-.2500*R4 

C(KONT)=0.00 

D(KONT)  =  -.5DO*TK4*RAY( I. (J*l ) . 1 ) ♦ ( . 5D0*T K4- . 2500 * R3 > *R A Y (  I, J.l) 

A  -»5*DS*(FLXB(J)*FLXL(I )) 

GO  TO  165 

213  CONTINUE 

:  CONVECT  ON  BOTH  SIOES 

IFLAG(4)=2 

IF( JFLG.NE.l)  GO  TO  2130 
E ACKONT ) =0 . DO 

EB(KONT)=-.5DO*TK4-.25DO*R4 

C(KONT)=.5DO*TK4 

0  (KONT) =-.500*TKl*RAY( (1-1) «J*1>*( .500*TK1*H( INDEX) *0S- .25D0*R3 > • 
A  RAY( I, J»1)-.5D0*H( INDEX > * ( TMPL ( I > *TNPB ( J > ) 

60  TO  165 
2130  CONTINUE 

EA(KONT)=.5DO*TK1 

EB(KONT)=-.5O0*TKl-.25D0*R4 
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92  2) 

923) 

924) 

925) 

926) 

927) 
928  ) 
929) 
930  ) 

931) 

932) 

933) 

934) 

935) 

936) 

937) 

938) 

939) 
94  0  ) 
94  1) 

942) 

943) 

944  ) 

945  I 
94b) 

947) 

948) 

949) 

950) 
951  ) 

952) 

953) 

954) 

955) 
95b) 
957) 

956) 

959) 

960) 
961  ) 

962) 
963  ) 

964) 

965) 
96b  > 

967) 

968) 

969) 

970) 

971) 

972) 

973) 

974) 

975) 

976) 

977) 

978) 

979) 
96  0) 

981) 

982) 

963) 

984) 

985) 

986) 
967) 

988) 

989) 

990) 

991) 

992) 

993) 

994) 

995) 

996) 

997) 

998) 

999) 
1000) 
1001) 
1002) 

1003) 

1004) 

1005) 

1006) 

1007) 

1008) 

1009) 

1010) 
1011) 
1012) 

1013) 

1014) 


C<KCNT>=C.D0 

D<KGNT)=-.5D0*TK4*RAY< I . ( J* 1 ) . 1 ) ♦ ( . 5 0 0 » T K4 *H (1 NDE X ) * DS- . 2 5D 0 *R 3  )  * 
&  R  A  Y  ( I»J»1)-.5D0*H( INDEX)* ITHPLC I)*THPB< J)  ) 

GO  TO  165 

214  CONTINUE 

SEMI-INF  ON  60TH  SIDES 
IFLAG(4>=3 
D I N  F  =  D 1  2/DI 

1 F  < JFLG.NE  .1 >  GO  TO  21 4C 
E  A  t  KONT >  =  0.D0 

E  H  C  KONT >=-TK4*DINF-CDI2*D12)*R4 
C (K0NT)=TK4»DINF 

0(K0NT)=-TRl*01NF.RAY«<I-l).J.l)-TKUtJ)*«DINF>* 

&  T  M  P  L  ( I )-TK  < I »  J  ) *GINF*  TMFb(J)* (TKi *OINF 
&  «2.DQ*TK<I«J>*DlNF-CDI2*DI2>*R3>*RAYCI,Jtl) 

GO  TO  165 
2140  CONTINUE 

E  A  t  KONT )=TK1*DINF 

EP  (KONT >  =  -TK 1*0  INF- CD1 2*012  >*R4 

C<KONT)=0.D0 

0(K0NT)=-TK( I, J) «  D INF* TMPL* I)-TK4*(DINF>* 

'&  RAY<I,CJ*1)»1)*<2.D0*TK1I,J)*CINF*TR4*D1NF- 
E  (DI2*OI2>*R3)*RAYU»J.l)-<TK(l,J)*DINF)*THPB(J) 

GO  TO  1 fc 5 

216  CONTINUE 

VERT  CONST  F  L  UX--H  OR  I Z  CONVECT 
1FLAG<4)=4 

IFC JFLG.NE.l)  GO  TO  216C 
CACKONT>=O.ODO 

EE<KONT)=-.5CD0*TK4-.250D0*R4 

C<KONT)=.bOOO*TK4 

O<KONT)=-.5OD0*TK 1*RAY  C  C 1-1 >. J,1)-.50D0*FLXL(I )*C$-.50D0* 

&  HUNUEX)*DS*TMPB(J)*(.50D3«TK3*.b0D0*HIINOEX)*DS-.250D0*R3)» 

&  R  A  Y  < 1«J«1) 

GO  TO  lb5 
2160  CONTINUE 

EAtKONT >=.50D0*TK  1 

Eh(KONT)=-.50DS«TKl-.5CD0*H(INOEX>*DS 

C<KONT)=0.CD0 

O(KONT>=-.50O0«TK4«KAY(l*tJ*l),l>-.50D0*FLXLCl >*DS-.50D0*HC INDEX) 
E  DS*TMPeCJ>*<.5CDO*TK4-.25300*R3>*RAYCI,J,I) 

GO  TO  165 

215  CONTINUE 

VERT  CONVECT--HOR 12  CONST  FLUX 
IFLAG(4)=7 

IFtJFLG.NE.l)  GO  TO  2150 
EA(KONT)=0.0D0 

EB«KONT)=-.5CD0*TK4-.50D0*HI I NDE X ) * DS- .250D0 *R 4 
C  (KONT)=.50D0*TK4 

OCKONT  >  =  -.2  5QOO*HClNDEX)*DS*TNPLCI ) -.5OD0 *FL XB I J ) *DS -.50DO *TK 1 
&  *RAY((I-1}»J.1)*(.50D0*TK1-«250D0*R3)*RAY<I»J»1) 

GO  TO  165 
215C  CONTINUE 

EAIKONT)=.50D0*TKl 

EBtKONT  >  =  -.50DC*TKl-.250D0*R4 

CCKONT>=0.000 

D(KONT)=-.50D0*TK4*RAY II »CJ-*1 > . 1 > - .50D0 *H C INDE X  )  *DS *TNPL ( I ) 

&  ♦C.50D3*TK4*.50DO*H(I  NDEX)*D$-.50DO*FLXB(J)*OS-.25000*R3) 

&  *RAY(I*J»1) 

GO  TO  165 

217  CONTI NUE 

VERT  CONST  FLUX--HOR I Z  SEHI-INF 
I  FLAG  ( 4 ) =5 

IFCJFLG.NE.l)  GO  TO  2170 
EACKONT )=O.DO 

EBIKONT )=-TK4*DI2-.2500*DI2*R4 
C(KONT)=TK4*DI2 

D (KONT >=-(TKl/(2. 00*01  >  >  *R  A  Y  <  <  I  - 1  >  ,  J  ,  1  >  -  ( TK  C  J  ,  J)  /  (2 . 00*D I  >  >  * 

&  TMPB(J)-FLXL<I)*DI2»DS*(TK1/(2.DO*OI)*TKII»J)/(2.DO*DI)- 
&  .25DO*OI2*R3)*RAY<I»Jtl> 

GO  TO  165 
2170  CONTINUE 

EAtKONT >=TKl/<2. 00*01) 

EB(KONT)=-TKl/(2. DQ*DI > -TK ( I » J > / ( 2 . DO *01 > - .250 0 »D 12 *R4 
C (KONT ) =0 • 00 

0(K0NT)=-TK4*0I2*RAY(I  «( J-»l)  » 1 )  -  ( TK  ( I  »  J)  /  <  2  .  DO  *01  >>*TNPB( J>- 
&  FLXLU )*DI2»DS*(TK4*OI2-.25DO*DI2*R3)«RAYCIfJ,l> 

GO  TO  165 

218  CONTINUE 

VERT=SEMI-INF»HORIZ=CONST  FLUX 
IFLAGC4)=9 

IFCJFLG.NE.l)  GO  TO  2180 
EAIKONT  >  =  0.000 

EB(K0NT)=-(TK4-*TMI*J)>/<2. 000*01)-. 5000  »D  12  *R  4 
C(KONT)=TK4/C2.0DO*OI ) 

D(KONT)=-FLXB(J)*OI2»DS-CTK(I»J)/(2.0DO*D1 ) ) *T NPL <1 > -TK 1 *01 2 
&  *RAY((I-l)*Jtl)**TKl*DI2-.50D0*DI2*R3)*RAY(I»J»l) 

GO  TO  165 
2180  CONTINUE 

EACKONT)=TK1*OI2 

EB( KONT >=-TK 1*012-. 5000*01 2*R4 

C (K  ONT )  =  0 . OD 0 

0CK0NT)=-CTK4/<2. 000*01 ) ) *R AY C I » I J*1 > » 1 ) -C TK Cl . J ) /< 2 . OD 0*0 I ) ) * 
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•  V ->  v 


1015)  &  THPL(I)-FLXB(J>*DI2*DS*((TK4*TK(ItJ))/(2.000*Dl)-.5000*OI2*R3> 

1016)  It  »R  A  Y  (  I  *  J  *  1 ) 

1017)  GO  TO  165 

1018)  219  CONTINUE 

1019)  C  VERT=CONVECTtHORIZ=SEMI-INF 

1020)  I  FLAG (  4 ) =6 

1021)  IF(JFLG.NE.l)  GO  TO  2190 

1022)  EA(KONT)=O.ODO 

1023)  EB(KONT)=-TK4*DI2-H(INDEX)*DI2-.50D0*DI2*R4 

1029)  C(K0NT)=TK4*DI2 

1025)  D(KONT)  =  -H( INDEX) *DI2*TMPL(  I ) - C TK 1 / C2 . OD 0*D 1 ) ) *RA Y! < 1-1 ) , J, 1 ) 

1026)  &  -(TK  (I  .J  ) /(2 .000  *DI )  >  *TMPB  (  J  >♦  (  (TK  1  +  TK  ( I  .  J)  )/ (2.  ODO  *01  > 

1027)  &  -.50DO*DI2*R3)*RAY(I.J.l) 

1028)  GO  TO  165 

1029)  2190  CONTINUE 

1030)  EA(K0NT)=TKl/(2. 000*01  ) 

1031)  £B(KONT)=-(TKl*TK(ItJ>)/(2.0D0*DI)-.50DC*DI2*R4 

1032)  C(KONT )=0.0D0 

1033)  D(KONT>=-H(INDEX>*DI2*TMPL(I)-(TK(I.J>/(2.0D0*DI>)*TMP8(J) 

1034)  &  -TK4*DI2*RAY(I,(J*1>. 1 ) ♦ < T K 4  *  0 1 2*H ( I  NOE  X ) *D 12 5 000 » DI 2* R 3 ) 

1035)  &  *R  AY  < 1 1  J» 1 ) 

1036)  GO  TO  165 

1037)  22 0  CONTINUE 

1038)  C  VERT  =S£MI-INF  t  HOR IZ=C0NVECT 

1039)  IFLAG(4>=8 

1040)  IFtJFLG.NE.il  GO  TO  2200 

1041)  EA(K0NT )=O.ODO 

1042)  EB(KONT)=-{TK4+TK(ItJ))/(2.0D0*Dl)-.50D0*DI2*R4 

1043)  C(KONT)=TK4/(2.0DC*DI) 

1044)  D(K0NT)=-H(IN0EX)*DI2*TMPB( J)-(TK(I.J)/(2.0DO*DI))*TMPL(I> 

1045)  &  -TK1*DI2*RAY((I-1).J.1)+(TK1*DI2*H(INDEX>*DI2-.50D0*DI2*R3) 

1046)  &  *RAY  (  1  *  J  *  1 ) 

1047)  GO  TO  165 

1048)  2200  CONTINUE 

1049)  EA(K0NT)=TK1»DI2 

1050  EB(KONT)=-TKl*DT2-H(INDEX)*DI2-.50D0*DI2»R4 

1051)  C (K  ONT ) =0 . ODO 

1052)  O(KONT>=-H(INOEX)*OI2*TMPB(J)-(TK4/(2.0Q0*OIl)*RAY(IttJ*l)tl) 

1053)  &  -(TK(ItJ>/(2.0O0*DI))*TMPL(I>'K(TK4*TK(ItJ>>/(2.0D0*DI> 

1054)  &  -.50D0*DI2«R3)*RAY(I,d,l) 

1055)  GO  TO  165 

1056)  211  CONTINUE 

1057)  C  CONST  CORNER  NODE 

1058)  EB(K0NT)=1 

1059)  0(K0NT)=RAY(I.J.1) 

1060)  IFLAG(4)=1C 

1061)  GO  TO  lb5 

1062)  133  CONTINUE 

1063)  C  BOTTOM  RIGHT-HAND  CORNER 

1064)  GO  TO  <2 21, 222. 223*224. 225.226*227, 228. 229*230). KRNRC3) 

1065)  222  CONTINUE 

1066)  C  CONST  FLUX  ON  BOTH  S10ES 

1067)  IF(JFLG.NE.l)  GO  TO  2220 

1068)  EA(K0NT )=.5D0*TK2 

1369)  EB(KONT)=-.5DO*TK2-.25DO*R4 

1070)  C(KONT)=O.DO 

1071)  D (KONT)=-.5O0*TKl *RAY( ( I -1 ) . J t 1 > ♦( . 5D0 *TK 1 - .25D0 *R3 > »R A Y ( I . J . 1 > 

1072)  &  -.500*DS»(FLXB(J)*FLXR(I)> 

1073)  GO  TO  165 

1074)  2220  CONTINUE 

1075)  EA(KONT)=.500*TK1 

1076)  EB(KONT)=-.5D0*TKl-.25DQ*R4 

1077)  C (KONT ) =0 . DO 

1078)  D(KONT)=-.5D0*TK2*RAY( i t (J-l)tl >♦< . 5D Q*TK2-. 25D0* R3 ) *R AY ( It J.l) 

1079)  &  -.5D0*US*(FLXB(J)*FLXR(I)> 

1080)  I FL  AG ( 6 ) =1 

1081)  GO  TO  165 

1082)  223  CONTINUE 

1083)  C  CONVECTIVE  ON  BOTH  SIDES 

1084)  IFLAG<6)=2 

1085)  IFIJFLG.NE.il  GO  TO  2230 

1086)  EA(KONT)=.5DO*TK2 

1087)  EB(KONT)=-.5D0*TK2-.25D0*R4 

1088)  C(KONT)=O.DO 

1089)  D (KONT )=-. 5D0*TK1*RAY( (1-1) » Jtl)*(.500« TK1*H( INDEX) ♦ OS- . 25D0*R 3 ) • 

1090)  &  RAY(ItJtl)-.5DO»H(INOEX)*(TMPB(J)4TMPR(I)) 

1091)  GO  TO  165 

1092)  2230  CONTINUE 

1093)  EA(KONT)=.500*TK1 

1094)  EB(KONT ) =- . 500 *TK 1 -.25 D 0 *R4 

1095)  D (KONT ) =0 .DO 

1096)  D(KONT)=-.5DO*TK2*RAY( It(J-l)*l)*(.5D0*TK2*H(INDEX) *DS-.2500*R3) » 

1097)  6  RAY(ItJtl)-.5D0*H(INDEX)*(TMPR(I)+TMPB(J)) 

1098)  GO  TO  165 

1099)  224  CONTINUE 

1100)  C  SEHI-INF  ON  BOTH  SIDES 

1101)  I  FLAG  ( 6 )  =  3 

1102)  0INF=0I2/DI 

1103)  IF(JFLG.NE.l)  GO  TO  2240 

1104)  EA(K0NT)=TK2*DINF 

1105)  EB(K0NT)=-TK2*DINF-(0I 2*012 )*R 4 

1106)  C(KONT)=O.DO 

1107)  D(K0NT)=-TK1*DINF*RAV( ( I -1 ) t J* 1 ) -TK ( I t J) *DI NF* 
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1108)  &  TMPR(1>-TK<I.J)«0INF*TMPB(J)*<TK1*DINF* 

1109)  &  2.DO*TK(I,J)*DINF-(OI2*DI2>»R3>*RAY(I.J»1) 

1110)  GO  TO  1 65 

1111)  2240  CONTINUE 

1112)  E  A  <  KONT )=TK1*DINF 

1113)  EB(K0NT)=-TK1»DINF-(DI2»DI2)*R4 

1114)  C ( KONT ) =C • DO 

111b)  D(K0NT)=-TK<I.J)*0INF«TMPR(I)-TK2«D1NF» 

1116)  &  RAY ( I .  <  J  - 1 )  « 1 >♦( 2.D0*TK (I *J)*DINF*TK2*DINF- 

1117)  &  DI2*D12»R3)*RAY< I »  J . 1 ) -  ( TK ( I*J)*DINF)*THPB(J) 

1116)  GO  TO  1 o5 

1119)  226  CONTINUE 

1120  C  VERT  CONST  FLUX--  HORIZ  CONVECT 


1121)  IFLAG<6):4 

1122)  IF<JFLG<NE<1 )  GO  TO  2260 

1123)  E  A  t  KONT ):.50D0*TK2 

1124)  EB (KONT):-. 5000 *TK2-< 250  DO *R4 

1125)  C<KONT)=0.0D0 

1126)  0<KONT)=-.50DO*TK1»RAY(<I-1)»J»1>-.5000*FLXR<I > *0 S- < 5000 »H( INDEX) 

112  7)  &  *DS*TKPB  < J) ♦<  <50DO*TK 1*.5CD0*H< INDEX) » DS- . 250 DO* R3) *RA Y< I , J, 1 ) 

1128)  GO  TO  165 

1129)  2260  CONTINUE 

1130)  E A(KONT):.50D0»TK 1 

1131)  EB<KONT>:-.50D0*TKl-.50DD*HCINDEX)*DS 

1132)  C  (KONT  ):0 .000 

1133)  D<KONT):-.50DO»TK2*RAY<1,CJ-1)«1)-.50DO*FLXRU )*DS-.50D0*H(INDEX) 

1134)  &  «DS*TMP8< J)*< .50O0*TK2-.250D0*R3> *RAY< I , J, 1 ) 

1135)  GO  TO  165 

1136)  225  CONTINUE 

1137)  C  VERT  CONVECT--  HORIZ  CONST  FLUX 

1138)  IFLAG(6)=7 

1139)  IF(JFLG.NE.l)  GO  TO  2250 

1140)  EA(KONT)=.50D0*TK2 

1141)  EB<KONT)=-.50DU*TK2-.50DC*H<INDEX)*DS-.250DC*R4 

1142)  C<KONT):0.0D0 

1143)  D(KONT):-.bOOO*H( INDEX  >*OS*THPR< 1 ) -.5  0 DO *FLXB < J ) *DS -<50D 0 *TK 1 

1144)  S  *RAY<<I-1).J»1)+(.50D0*TK1-.25DD0*R3>*RAY(I»J»1) 

1145)  GO  TO  lbb 

1146)  2250  CONTINUE 

1147)  EA<KONT):.50D0*TKl 

1148)  EB(KONT )=-.5000*TKl-.250DO*R4 

1149)  C <KONT)=O.ODO 

1150)  D(KONT)=-.5CDO*TK2*RAY (I «  t  J-l ) *1 )-.50DO»HCINDEX)*DS*THPR(I )♦ 

1151 )  &  < .5000*TK2*.5000*H< I N OE X ) *DS- .50D0 *FLXB < J) *DS-.250D0*R3) 

1152)  &  »RAY  < I  *  Ji 1 ) 

1153)  GO  TO  165 

1154)  227  CONTINUE 

1155)  C  VERT  CONST  FLUX*  HORIZ  SEHI-INF 

1156)  IFLAG<6)=5 

1157)  IF(JFLG.NE.l)  GO  TO  2270 

1158)  EA(K0NT)=TK2*DI2 

1159)  EBtKONT ):-TK2*DI2-.25D0*DI2*R4 

1160)  C<KONT):O.DO 

1161)  D(KONT)=-(TKl/(2.O0*0I))*RAY((I-l)«Jtl)-(TK(I.J)/(2«D0*DI))* 

1162)  &  TMPB<J)-FLXR(I)*OI2*DS*(TKl/t2.DO»DI )»TK(I»J)/(2.D0*DI >- 

1163)  &  .25D0*DI2*R3)*RAY(I.J»1) 

1164)  GO  TO  165 

1165)  2270  CONTINUE 

1166)  EA(KONT)=TK1/(2.DO*DI) 

1167)  EB<KONT ):-TKl/<2. 00*01 )-TK( I»J)/(2.D0*DI )-.25D0*DI2*R4 

1168)  C(KQNT)=0.00 

1169)  0<K0NT):-TK2*DI2*RAY(I . < J-l ) . 1 ) - < T K < I « J ) / < 2 . DO *D I ))*TMPB(J>- 

1170)  &  FLXR(I)»OI2*OS-»<TK2»D12-.25DO*DI2*R3)*RAY(I»J»1) 

1171)  GO  TO  165 

1172)  228  CONTINUE 

1173)  C  VERT=SEHI-INF.HORIZ=CONST  FLUX 

1174)  I  FLAG  <  6 ) :9 

1175)  IF(JFLG.NE.l)  GO  TO  2280 

1176)  EA<KONT):TK2/(2.0U0»DI  ) 

1177)  EBtKONT ):-(TK2*TK  < 1 » J ) ) /(2.0DO*DI)-.50D0*D12*R4 

1178)  C(KONT)=O.ODO 

1179)  D<KONT)=-FLXB(J>*DI2*DS-<TK(I*J)/C2.0D0*DI>)*THPRCI) 

1180)  &  -TK1*DI2*RAY<(I-1)»J»1)*<TK1»DI2-.50D0*DI2*R3)*RAY(I»J»1) 

1181)  GO  TO  165 

1182)  2280  CONTINUE 

1183)  EA(K0NT)=TK1*0I2 

1184)  EB(KONT):-TK1*DI2-.50DO*DI2*R4 

1165)  C(KONT)=O.ODO 

1186)  0<KONT)=-<TK2/(2.0DO*DI))*RA Y < I , < J-l ) , 1 ) -< TK < I . J) / t 2 . ODO* DI ) ) * 

118  7)  &  TMPR  < I )-FLXB< J)*DI2*DS* < <TK2*TK  <1 . J) )/<2.0D0*DI )-.50D0*DI2*R3)» 

1188)  &  R A Y ( I,J,1) 

1189)  GO  TO  165 

1190)  229  CONTINUE 

1191)  C  VERT=CONVECT»HORIZ=SEHI-INF 

1192)  I  FLAG ( 6 ) =6 

1193)  IF(JFLG.NE.l)  GO  TO  2290 

1194)  EA(K0NT)=TK2»0I2 

1195)  EB < KONT )=-TK2«DI2-H< INDEX) •DI2-.50D0*DI2*R4 

1196)  C(KONT>=0.000 

1197)  0(K0NT)=-H(lNDEX)*DI2*TMPR(I)-(TKl/(2.000*DI))*RAY((I-l)tJ»l) 

1198)  (  -(TK (I «J )/<2.0D0*DI ) ) *THPB(J)*( CTKl+TKtltJ) )/ <2.0D0*DI ) 

1199)  C  -.5000*DI2»R3)»RAY(IfJ»l) 

1200)  GO  TO  165 
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• 

m 

1201) 

2290 

CONTINUE 

1202) 

EAIKONT)=TK1/I2.0DO*DI) 

1203) 

EB IKONT >=-ITKl*TK  II  ,J> ) / 12 . 000 «D I ) - .5 OOO *0 I 2 *R 4 

1204) 

C(KONT)=O.ODO 

1205) 

OIKONT)=-HI INDEX)*0I2«THPRI I ) - 1 T K 1 1 . J ) / 1 2 . 000* 0 1 >  ) *THPB I J ) 

.  * 

1206) 

l  -TK2*OI2*RAYII,IJ-1),1)4<TK4*OI24HIINDEX)*DI2-.50DO*DI2*R3)* 

*  * 

1207) 

l  RAYII.J.1) 

1206) 

GO  TO  165 

a 

1209) 

230 

CONTINUE 

a 

1210) 

C 

VERT=SEMI -INF,H0RI2=C0NVECT 

1211) 

I  FLAG ( 6 ) =8 

,  - 

1212) 

IFIJFLG.NE.l)  GO  TO  2300 

1213) 

EA(KONT)=TK2/<2.0DO*D1  ) 

•»  * 

1214) 

EB(K0NT>=-(TK2*TK (I.J) ) / ( 2 . D0*DI ) - .50D0 *01 2 *R4 

1215) 

C IKONT )  =  0  «  ODO 

- 

1216) 

DIKONT)=-HI INDEX) *0 12* THPB 1 J > -  1 TK(1 ,J)/I 2. 000*01 ) )*TNPRII ) 

1217) 

&  -TK1*D12*RAY II I-1),J,1)4(HI INDEX) *DI24TK1*OI2-.50DO*DI2*R3)* 

r\ 

1218) 

&  RAYlItJtl) 

a 

1219) 

GO  TO  165 

■ 

1220) 

2300 

CONTINUE 

1221 ) 

EAIK0NT)=TK1*0I2 

1222) 

EBIKONT>=-TK1*DI2-.50DO*DI2»R4-H{  INDEX) *012 

1223) 

C ( KONT) =0 • 000 

k 

1224) 

OlKONT)=-HIINDEX)*DI2*THPB»J»-iTK2/(2.0DO*DI))*RAYtl.l J-l),l> 

1225) 

&  -ITKII,J)/I2.0D0*DI))*THPRII)4| 1 T K 24TK 1 I ♦ J ) > / 12 . ODO *DI ) 

Sv 

1226) 

&  -.5000*DI2*R3>*RAYII,J.l) 

1227) 

GO  TO  165 

1228) 

221 

CONTINUE 

m 

1229) 

C 

CONSTANT  CORNER  NODE 

1230) 

EB( KONT ) =1 

1231 ) 

DiKONT)=RAYIItJtl) 

V 

1232) 

I  FLAG  I6)  =  l 0 

1233) 

GO  TO  165 

./ 

1234) 

104 

CONTINUE 

s,< 

1235) 

C 

TOP  RIGHT-HANO  CORNER 

'< 

1236) 

GO  TO  1231, 232, 233, 234, 235,236, 237, 238,239, 240), KRNRI4) 

.*' 

1237) 

232 

CONTINUE 

. ' 

1238) 

C 

CONST  FLUX  ON  BOTH  SIDES 

Ss 

1239) 

IFIJFLG.NE.l)  GO  TO  2320 

N 

1240) 

EAIKONT)=.5D0*TK2 

P 

1241  > 

EBIKONT)=-.5DO*TK2-.25DO*R4 

m  * 

1242) 

CIKONT)=O.DO 

‘.  * 

1243) 

D IKONT )  =  -.5D0*TK3 *R AY  1 II *1 ) « J, 1 ) ♦( . 5D 0 *T K3-. 25D0*R3 ) * 

• 

1244) 

&  RAYl I,J,1)-.5D0*DS*IFLXTI J)4FLXRII )) 

i  *t 

1245) 

GO  TO  165 

*  \ 

1246) 

2320 

CONTINUE 

*  *. 

1247) 

EAIKONT)=0.D0 

V*. 

1248) 

EBIKONT)=-.500*TK3-.250G*R4 

al 

1249) 

C<KONT)=.500*TK3 

■ 

1250) 

DIKONT)=-.5O0*TK2*RAY< I,IJ-1),1 ) ♦  ! . 5D0*TK2 -,25D0*R3) *RA Y 1 I , J ,1 ) - 

m 

1251) 

4  .5D0*DS*IFLXT IJ)*FLXR II )) 

1252) 

IFLAGI8)=1 

v 

1253) 

GO  TO  165 

1254) 

233 

CONTINUE 

1255) 

C 

CONVECTIVE  ON  BOTH  SIDES 

. ", 

1256) 

IFLAGI8)=2 

,N 

1257) 

IFIJFLG.NE.l)  GO  TO  2330 

r* 

1258) 

EAIKONT)=.5D0*TK2 

1259) 

EBIKONT)=-.5D0*TK2-.25D0*R4 

■ 

1260) 

CUONT)  =  0.D0 

1261) 

0<K0NT)  =  -«500*TK3*RAY((1«1)»J»1)*(.5D0*TK3*HII NOEX) »  OS- .2500*R3 ) * 

*  - 

1262) 

l  RAYII,J»1>-.5DO*HIINOEX)*TNPT(J)-.5DO*HIINDEX)*THPR II) 

L  . 

1263) 

GO  TO  165 

•  • 

1264) 

2330 

CONTINUE 

l.'- 

1265) 

E A IKONT )=Q .DO 

it* 

1266) 

EBI KONT )=-. 500 *TK 3-.25D0 *R4 

,.•' 

1267) 

CIKONT)=.5DO*TK3 

>  ' 

1268) 

OIKONT )  =  -.500*TK2*R A Yl  It IJ-1),1 ) 41.50 0*TK24H( INDEX) *DS-.2500*R3)* 

*4 

1269) 

4  RAVII,J,1)-.5D0*HIINDEX)*TMPTIJ)-.5D0*HIINDEX)*TMPRII) 

y 

1270) 

GO  TO  165 

P 

1271) 

234 

CONTINUE 

1272) 

C 

SEMI-INF  ON  BOTH  SIOES 

►  * 

1273) 

I  FLAG  1 8) =3 

<> 

1274) 

WRITE! 1 *9997 )  I.J 

1275) 

GO  TO  9999 

r  * 

1276) 

236 

CONTINUE 

1277) 

C 

VERT  CONST  FLUX--HOR I Z  CONVECT 

1278) 

I  FLAG  18)  =  4 

1279) 

IFIJFLG.NE.l)  GO  TO  2360 

1280) 

EAIKONT)=.5000»TK2 

R 

1281 ) 

EB IKONT )=-.50D0«TK2-.250O0»R4 

1282) 

CIKONT)=O.ODO 

1283) 

OIKONT)  =  -.50O04TK3«RAY II l4l),J,l)-.50DO*FLXRII  )*DS 

1284) 

6  -.50D0*HI INDEX )»DS*TNPTIJ) 4 |.50D0»TK3*.50D0*H< INDEX )*DS 

».i 

1285) 

l  -.25000*R3)*RAYII,J,1) 

».1 

1286) 

GO  TO  165 

.  7 

1287) 

2360 

CONTINUE 

%* 

1288) 

EAIKONT)=O.ODO 

i \ 

1289) 

EBIKONT >=~. 5000 *TK3-.5000*HI INDEX )*OS 

A 

1290) 

CIKONT)=.SOOO*TK3 

m 

1291) 

OIKONT)=-.5000*TK2*RAYIItlJ-l)tl)-.5000*FLXRII )*DS 

»/• 

129  21 

1  9Q  *  1 

l  -.5000«HI INDEX )*DS*THPTIJ)4«.50D0*TK2-.250D0*R3)»RA YII, J,l) 

a  * 

1294) 

C 

VERT  SIOE=SEMI-INF,  HORIZ  CONST  FLUX 

*.* 

50 

“i 

i 

'.■  V  W  v 
,>*.%>  />  . 

a 

:  v  %v.  v-v 

„V%V,_v 

V  * 
^ 

1295) 

1296) 

1297) 

1298) 

1299) 
1300  ) 

1301 ) 

1302) 

1303) 

1304) 

1305) 

1306) 

1307) 

1308) 

1309) 

1310) 
1311  > 

1312) 

1313) 

1314) 

1315) 

1316) 

1317) 
1318  ) 

1319) 

1320) 

1321) 

1322) 

1323) 

1324) 

1325) 

1326) 

1327) 

1328) 

1329) 

1330) 
1331  ) 

1332) 

1333) 

1334) 

1335) 

1336) 

1337) 

1338) 

1339) 

1340) 

1341) 

1342) 

1343) 

1344) 

1345) 

1346) 

1347) 

1348) 

1349) 

1350) 

1351) 

1352) 

1353) 

1354) 

1355) 

1356) 

1357) 

1358) 

1359) 

1360) 

1361 ) 

1362) 

1363) 

1364) 

1365) 

1366) 

1367) 

1368) 

1369) 

1370) 

1371) 

1372) 

1373) 

1374) 

1375) 

1376) 

1377) 

1378) 

1379) 

1380) 

1381) 

1382) 

1383) 

1384) 

1385) 

1386) 

1387) 

1388) 


238  CONTINUE 
IF(JFLG.NE.l)  GO  TO  2380 
EA(KONT )=TK2/(2.0D0*DI  ) 

EB(KONT)=-.5000*(TK(ItJ)/DI*TK2/DI*DI2*R4> 

C(KONT)=O.ODO 

D (K0NT)=-TK3*D12*RAY( ( I ♦ 1 > » J »1 ) -FI  XT < J > *DI 2 *DS ♦ 

&  <TK3»DI2-.50O0*DI2*R3)*RAYI If Jt 1) 

GO  TO  165 
2380  CONTINUE 

EA(KONT)=O.ODO 

Ee<KONT>=-TK3»DI2-.50DO*OI2*R4 

C(KONT)=TK3«0I2 

D (KONT)-(-TK(I , J) /(2.0D0*DI ) >»TMPRU  >-<TK2/«2.0D0*Dl ) )• 

&  RAY ( If ( J  +  l ) f 1 ) ♦ <  T  K  < 1 , J) /<2.0DO*DI)*TK2/(2.0DO*DI >-.50D0*DI2 
It  »R3)  »R  AY  (  I  fjfl  > - FLXT  <J)*DJ2*DS 
GO  TO  165 
235  CONTINUE 

C  VERT  CONVECT --HOR 1 2  CONST  FLUX 

IFLAG(e>=7 

IFtJFLG.NE.l)  GO  TO  2350 
EAIKONT )=.50DO*TK2 

EBIKONT  >=-.50D0»TK2-.50D0*H( I NDEX > *DS - .25000 *R 4 
C (KONT  >  =0 . 0D0 

DIKONT)=-.50D0*FLXT(J)*OS-.5  0D0»T<3«RAY(U*l)fJfl>* 

&  (.50D0*TK3-.250D0*R3>*RAY(IfJtl>-.5QD0*H(INDEX)»DS*TMPMI> 

GO  TO  165 
2350  CONTINUE 

E  A (KONT )  =  0 • 0 DO 

E6(KGNT)=-.50O0*TK3-.25GD0*R4 
C (K0NT)=.50DG*TK3 

D(KONT)=-.50D0*TK2»RAY ( I « < U-l ) , 1 ) - . 50 DO *H C I NDE X > »DS* 

&  TMPR (I )♦ ( .5  000 *TK2*.5  ODO*H( INDEX > *DS-. 5 CDO *FL XT ( J> *DS 
&  -.25COO»R3)*RAY(IfJfl) 

GO  TO  165 
240  CONTINUE 

C  VEKT=SEMI-INF,HORIZ=CONVECT 

IFLAG(8)=8 

IFtJFLG.NE.l)  GO  TO  240C 
E  A ( KONT )  =  Tk2/(2.0DC*D1  ) 

EB(KONT)=-(TK2+TK<IfJ))/(2.0D0»DI)-.50D0*DI2*R4 
C ( KONT )  =  0 . CDO 

D(KONT)=-H(INDEX) «012*TKPT4J)-CTK(I fJ)/C2.000*Dl) >*TMPR(I> 

&  -TK3*DI2»RAY4(I*l)fJfl)*(H(INDEX)»DI2*TK3»DI2-.50D0*DI2*R3>* 

&  RAY(lfJtl) 

GO  TO  165 
2400  CONTINUE 

EACKONT )=0.0D0 

EB(KONT)=-TK3*DI2-.50D0*DI2«R4-H(INDEX)*DI2 

CCK0NT)=TK3*DI2 

0 (KONT )=-H (INDEX) «0I2*TMPT (J)-(TK2/(2.0DC*DI >)*RAY(I ,CJ-1 >,1) 

&  -(TK(IfJ)/(2.0D0*DI))*TKPR(I)*((TK2>TK(IfJ))/(2.0D0*DI)-.50D0* 
&  0I2*R3)*RAY(IfJfl) 

GO  TO  165 
237  CONTINUE 

239  CONTINUE 

WRI TE ( 1 t 9997 )  I,J 
GO  TO  9999 
231  CONTINUE 

C  CONSTANT  CORNER  NODE 

EBCKONT)=l 
D(KONT)=RAY(If J»l) 

IFLAG(8)=1C 
GO  TO  165 

C  ((((((((((((((((((((((((((((I 


165  CONTINUE 

C  CHECK  WHICH  PASS 

IF< JFLG.E0.2)  GO  TO  1602 
170  CONTINUE 

C  A  ROW  SET  UP  ON  1ST  PASS  NOW. 

L  =  X 
I  F=  1 

CALL  TRIOIGtLtEAf EBtCf D) 

DO  1607  K=ltX 
T  EMP  C I f K ) =D ( K  > 

1607  CONTINUE 
160  CONTINUE 

C  THE  WHOLE  GRID  1ST  PASS  COMPUTED  NOW. 

JFLG=2 
KONT=0 

DO  173  I=lfY 
00  174  J=lfX 
OLDTIIt J)=RAY(I, J,l) 

RAY(IfJfl) =  T  EMP ( I f J) 

174  CONTINUE 
173  CONTINUE 

GO  TO  2001 
1602  CONTINUE 
1700  CONTINUE 

C  A  COLUMN  SET  UP  ON  2ND  PASS  NOW. 

1613  CONTINUE 
IF=1 
L  =  Y 
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1%70>  404  FORHAT(/*2X, ‘BOUNDARY  CONDITIONS:*) 

1471)  DO  410  K  =  l«8 

1472 )  KKK=K 

1473)  IF(IFLAG(K).EQ.O>  GO  TO  410 

1474)  GO  TO  (400*450»401»451»402*452*403*453)*KKK 

1475)  400  WRITE (6*405) 

1476)  405  FORMAT (/»1X» ‘TOP  BOUNDARY*) 

1477)  60  TO  470 

1478)  401  WRITE (6*406) 

1479)  406  FORMAT (1X»*LEFT  BOUNDARY*) 

1480)  60  TO  470 

1481)  402  WRITE(6*407) 

1482)  407  FORMAT (IX* ’BOTTOM  BOUNDARY*) 


v.'.  ■  ■ 


.  .  .s . 
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1483) 

1484) 

403 

1485) 

408 

1486) 

1487) 

450 

1488) 

460 

1489) 

1490) 

451 

1491) 

461 

1492) 

1493) 

452 

1494) 

462 

1495) 

1496) 

453 

1497) 

463 

1498) 

1499) 

470 

1500) 

1501) 

1502) 

471 

1503) 

492 

1504) 

1505) 

472 

1506) 

493 

1507) 

1508) 

473 

1509) 

494 

1510) 

1511) 

474 

1512) 

484 

1513) 

1514) 

475 

1515) 

485 

1516) 

1517) 

476 

1518) 

486 

1519) 

1520) 

477 

1521) 

487 

1522) 

1523) 

478 

1524) 

488 

1525) 

1526) 

4  79 

1527) 

489 

1528) 

1529) 

480 

1530) 

491 

1531) 

1532) 

481 

1533) 

1534) 

1535) 

482 

1536) 

1537) 

1538) 

483 

1539) 

1540) 

1541) 

416 

1542) 

417 

1543) 

418 

1544) 

410 

1545) 

1546) 

1547) 

183 

1548) 

1549) 

1550) 

1551) 

1552) 

C 

1553) 

C 

1554) 

1555) 

1556) 

511 

1557) 

1558) 

512 

1559) 

1560) 

513 

1561) 

1562) 

514 

1563) 

1564) 

515 

1565) 

i567) 

516 

1568) 

1569) 

\l7A\ 

519 

1572) 

1573) 

520 

157$) 

521 

1576) 

.'* 

V-V-* 

60  TO  470 
WRITE (6*406) 

F0RHAT(1X«*RI6HT  BOUNDARY*) 

60  TO  470 
WRITE  (6*460 ) 

FORMA TC IX » • TOP  LEFT  CORNER*) 

GO  TO  470 
WRITE (6*461 ) 

F OR HATC1X** BOTTOM  LEFT  CORNER*) 

GO  TO  470 
W  R ITE ( 6*  462 ) 

F OR MAT ( IX* * BOTTOM  RIGHT  CORNER*) 

GO  TO  470 
WRITE (6*463) 

FOR  MAT ( IX  *  *  TOP  RIGHT  CORNER*) 

GO  TO  470 
CONTINUE 

GO  TO  (471 *472 *47 3 *474 *475 *4 76 *4 77 *478 *4 79 *480 *48 1*482* 
6  483) « IFLAG(K ) 


WRITE (6*492) 

FORMAT (1H>«*  CONST  FL 

GO  TO  410 
W  R I TE  (6  *493 ) 

FORMAT! 1H*« •  CONVECTIVE 

GO  TO  410 
WRIT£(6»494) 

F  ORMAT  (  1H*-*  •  SEMI-INF 

CO  TO  410 
wr*  I  TE  (6*484) 

F Of* MAT ( 1H*. •  VERT  CON 

GO  TO  410 
WRITE (6*485) 

FORMAT ( 1H** •  VERT  CONST 

GO  TO  410 
WRITE (6*486) 

FORMAT  (1H-**  *  VERT  CO 

GO  TO  410 
WRITE(6*487) 

FORMAT ( 1H-*,  •  VERT  CON 

GO  TO  410 
WRI TE (6*488  > 

FORMAT ( 1H*» *  VERT  SEMI 

GO  TO  410 
WRITE (6*489) 

F ORMAT ( 1H-*.  •  VERT  SEMI- 

GO  TO  410 
WRITE(6*491) 

FORMAT ( 1H+* •  CONST) 

GO  TO  410 
CONTINUE 
WRITE (6*416) 

GO  TO  410 
CONTINUE 
WRITE (6*417) 

GO  TO  <10 
CONTINUE 
WRITE (6*418) 

GO  TO  410 

FORMAT ( 1H*, •  SEMI-INFINITE*) 

FORMAT ( 1H*« •  CONVECTIVE  SURFACE*) 

FORMAT ( 1H*« •  CONSTANT  HEAT  FLUX*) 

CONTINUE 
WRITE (1*183)  Q 
W  R I TE (6  *  183 )  Q 

FORMAT (/» IX* ‘NUMBER  OF  TIME  STEPS=»»I5) 

00=0*DELT 
RRX=X/17 
JJXsRRX 

JJR=X-(17*JJX) 

WRITE  THE  FINAL  TEMPS  t  THERMAL  VALUES 
00  510  KK=1 *5 

GO  TO  (511*512*513, 514*515) *KK 
WRI TE (6*893 )  Q,DD 
GO  TO  516 
WRITE (6*911) 

GO  TO  516. 

WRITE (6*912) 

GO  TO  516 
WRITE (6*914) 

GO  TO  516 
WRITE (6*913) 

GO  TO  516 
DO  518  JJsl*UJX 
IF(JJX.EQ.O)  GO  TO  518 
J2=17*JJ 
Jl= J2-16 

GO  TO  (519*520*521*522*523) «KK 
WRITE (6*500)  ((RAY(I*J,1),J:J1,U2>*I:1,V) 

GO  TO  517 

WRITE (6,500)  ((TK(I*J) *J=J1*J2)*I=1«Y) 

GO  TO  517 

WRITE (6 .500)  ((CP(I«J),J=Jl*J2)*Isl*Y) 


CONST  FLUX  ON  BOTH  SIDES*) 


CONVECTIVE  ON  BOTH  SIDES*) 


SEMI-INF  ON  BOTH  SIDES*) 


VERT  CONST  FLUX— HORI2  CONVECT*) 


VERT  CONST  FLUX—  HORIZ  SEMI-INF*) 


VERT  CONVECT— HORIZ  SEMI-INF*) 


VERT  CONVECT— HORIZ  CONST  FLUX*) 


VERT  SEMI-INF— HORIZ  CONVECT*) 


VERT  SEMI-INF--HORIZ  CONST  FLUX*) 


CONSTANT  TEMPERATURE*) 


WKl It  «6*3( 

60  TO  517 


-  *  *  ’  “  m  m  *  ’  •  v  '  .  "  -  '  - 

■  •  V  V  *,*  v  •  .  -f. 


1577) 

1578) 

1579) 

1580) 

1581) 

1582) 

1583) 

1584) 

1585) 

1586) 

1587) 

1588) 

1589) 

1590) 

1591) 

1592) 

1593) 

1594) 

1595) 

1596) 

1597) 

1598) 

1599) 

1600) 
1601) 
16C2  ) 

1603) 

1604) 

1605) 

1606) 

1607) 

1608) 

1609) 

1610) 
1611) 
1612) 

1613) 

1614) 

1615) 

1616) 

1617) 

1618) 

1619) 

1620) 
1621) 
1622) 

1623) 

1624) 

1625) 

1626) 

1627) 

1628) 

1629) 

1630) 

1631) 

1632) 

1633) 

1634) 

1635) 

1636) 

1637) 

1638) 

1639) 

1640) 

1641) 

1642) 

1643) 

1644) 

1645) 

1646) 

1647) 

1648) 

1649) 

1650) 

1651) 

1652) 

1653) 

1654) 

1655) 

1656) 

1657) 

1658) 

1659) 

1660) 
1661  > 
1662) 

1663) 

1664) 

1665) 

1666) 

1667) 

1668) 
1669) 


WR I TE (6,53  G )  (IROUiJ) 
GO  TO  517 

WRITE(6,502)  ((1ST  AT (I 
IF(JJ.LT.JJX)  WRITE (6« 
CONTINUE 
J2  =  X 

J1=X-JJR*1 

IF(JJR.EQ.O)  GO  TO  524 
IF(JJX.NE.O)  UK1TE(6*5 
DO  524  1=1, Y 
GO  TO  (525,526,527,528 
WRITE(6,500>  (R AY ( I « J, 
GO  TO  524 

WRITE (6,500)  (TK(I«J)« 
GO  TO  524 

WRITE(6«500)  (CP(l,d), 
GO  TO  524 

WRITE(6,500)  ( R 0 (  I  ,  J)  , 
GO  TO  524 

WRITE  (6,502)  (1ST  AT (I, 

CONTINUE 

CONTINUE 

F0RHAT(1X,17F7.2) 
FORMAT!/, IX, ’FINAL  TK  ( 
FORMAT!/, IX, ’FINAL  CP( 
FORMAT! 1X,17F7.2) 
FORMAT!/, IX, ’FINAL  R0( 
FORMAT!/, IX, ’FINAL  1ST 
FORMAT (IX, 171 3) 
FORMAT!/, 3X, ’AND  THE 
WRITE(11,915)  CPPC(l), 
WR I TE ( 1 1 *916 )  TPC  (1),T 
F0RMAT(/,1X,’CPPC(1)=* 
F0KMAT(1X,’TPC(1)=’,F6 
CONTINUE 
FORMAT (/,/) 

CREATE  A  NEW  DATA 
WRITE(9,10)  DS ,OELT ,DI 
WRI TE (9,20 )  A, X , Y  *NIS0 
WRITE(9,16>  ((ROPC(K), 
WRITE(9,15) (TISO(B)»B= 
WRITE (9 ,30 )  (H(L),L=1, 


(  ( RO ( I , J) ,J  =  J1,J2),I  =  1,Y) 


,J),J=J1,J2),I=1,Y) 

503) 


, 529  > , KK 
1) ,J=J1, J2) 


(TK(  I «J>« J  =  J1,U2) 
(CP( I,J)» J=J1,J2> 

(  R  0 ( I,J),J=J1«J2) 

(  1ST  ATd,  J)  ,J=J1,  J2> 


i,j>:*> 

i,j):») 

l, j):*> 

AT(i,j):») 


NEXT  SET  OF  COLUMNS:’) 
ROPC(l) 

DEL 

, F6.2 , •  R0PC(1)=’,F6*2> 

.2,’  TDEL=*,F6»2> 


CREATE  A  NEW  DATA  FILE  WITH  END  RESULTS 
WRITE(9,10)  DS ,OELT ,DI »TDEL 
WRITE (9,20)  A,X,Y,NISO,ITRT ,IMAX,ITPC 
WRITE (9, 16)  ( (ROPC(K) ,CPPC(K),HL(K) , TPC ( K ) ) ,K=1,A) 
WRITE(9«15) (TIS0(6)*B=1«N1S0) 

WRITE (9 ,30 )  (H (L ) *L  =  1 «  A ) 

WRITE (9,35)  ((RAY(I*J«3)«J=1*X),I=1*Y) 

WRITE (9, 35)  ((RAY(I,J,2),J=1,X),I=1»Y) 

WRITE (9, 35)  ((RAY(I,J,1),J=1,X),I=1,Y) 

WRITE(9,34)  (KRNR(J)«J=1,4) 

WRITE(9,35)  (FLXT ( J),J=1 ,X) 

WRITE (9, 35)  (FLXB(J), J=1,X) 

WRITE (5,35)  (FLXL (I),I=1,Y> 

WRITE(9,35)  (FLXR(I ),I=1,Y) 

WRITE (9, 35)  (TMPT (J),J=1,X) 

WRITE (9,35)  (TMPB(J), J=1,X) 

WRITE (9,35)  (TMPL(1),I=1*Y) 

WRITE(9,35)  ( TMPR (I),I=1,Y) 

FORMAT ( IX ,’ NO  EQUATION  FOR  R AY ( » , 1 2 , » , • , 12, • ) . * ) 
CONTINUE 

CLOSE  FILES 

CALL  CONTRL(4,’ADPTMP’,ll) 

CALL  C0NTRL(4,’A0PDAT’,5) 

CALL  C0NTRL(4»’ADP0UT»»6) 

CALL  C0NTRL(4,’ADPNDT* ,9) 

CALL  C0NTRL(4»»P0INT1’,13) 

CALL  EXIT 


C  _ 

c  *  *  * 

C  TH 


SUBROUTINE  TR 1 0 1 G ( N, A , B, C*D> 

IMPLICIT  INTEGER* 2 ( I »N ) 

IMPLICIT  DOUBLE  PR E Cl S ION ( A *B,C ,D) 

DIMENSION  A(1),B(1),C(1),D(1) 

ADJUST  COEFFS  FOR  BAD  FROM  ELIMINATING 
DO  10  1=2, N 
AB=  A  (  I ) /Bd -1 ) 

B(I)=B(I)-AB*C(I-1) 

D(I)=D(I)-AB*D(I-1) 

CONTINUE 

BACK  SUBSTITUTE 

N1=N-1 

D(N)=D(N)/B(N) 

00  20  1=1, N1 

N  *  I 

D(M)=(D(M)-C(M)*D(M*1) )/B(M) 

CONTINUE 

RETURN 

END 


7TTTTT****** ************************* - 

SUBROUTINE  ISOTHM 

IS  PROGRAM  FINOS  ISOTHERMS  IN  RAY(ItJ*l). 
IMPLICIT  INTEGER*2(A,g«I-L«N,a*U*X*Y,Z) 
IMPLICIT  DOUBLE  PRECIS I0N€C-H,N*0*P*R-V > 
DIMENSION  EXRY(9*200) *  EYRY(9,200) ,COUNT (9) 
C0MM0N/M12I  /  X, Y,NIS0,0 
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1670) 

1671) 

1672) 

1673) 
1674  > 
167b) 

1676) 

1677) 
167S) 

1679) 

1680) 
1681  > 
1682) 

1683) 

1684) 

1685) 

1686) 

1687) 

1688) 

1689) 

1690) 

1691) 

1692) 

1693) 

1694) 

1695) 

1696) 

1697) 

1698) 

1699) 

1700) 

1701) 

1702) 

1703) 

1704) 

1705) 

1706) 

1707) 

1708) 

1709) 

1710) 

1711) 

1712) 

1713) 

1714) 

1715) 

1716) 

1717) 

1718) 

1719) 

1720) 

1721) 

1722) 

1723) 

1724) 

1725) 

1726) 

1727) 

1728) 

1729) 

1730) 

1731) 

1732) 

1733) 

1734) 

1735) 

1736) 

1737) 

1738) 

1739) 

1740) 

1741) 

1742) 

1743) 

1744) 

1745) 

1746) 

1747) 

1748) 

1749) 

1750) 

1751) 

1752) 

1753) 

1754) 

1755) 

1756) 

1757) 

1758) 

1759) 

1760) 


C0MM0N/H12R /  R A Y < 80 ,80 , 3 ) ,T I  SO ( 9 ) ,DS, DELT 
C  NISO=NUMBER  OF  ISOTHERMS 

C  COUNT  <B)=COUNTER  FOR  NO.  ELTS  IN  EACH  ISOTHERM 
C  T ISO(B) =  TEMP  OF  EACH  ISOTHERM 

C  _  _  LET  TISO(l)  BE  THE  HOTTEST  ISOTHERM 
C  EXRYIB»K)=  ARRAY  FOR  X-C OORD I  NATES 
C  B  INDICATES  WHICH  ISOTHERM 

C  KrCOUNT ( B)  l  INDICATES  POSITION  OF  PT  IN  ISOTHERM  LIST 

C 

C 

C  CHANGE  THE  NEXT  STATEMENT  TO  AGREE  WITH  DIMENSION 
DO  6  K=l,200 
DO  7  B  =  1 «  N I  SO 
E  XR  Y ( B t  K ) =0 
EYRY(B*K)=G 

7  CONTINUE 

6  CONTINUE 

DO  8  K=1 »9 
COUNT (K )=0. 

8  CONTINUE 
CNTPT=0. 

C  LOOP  TO  LOCATE  ISOTHERMS 
C  FOR  LEFT  SEMI -INNF  .LET  XL  =  3 

C  FOR  BOTTOM  SEMI-INF, LET  YB=Y-3 

C  FOR  RIGHT  SEMI-INF, LET  XR  =  X-3,  AND  COMMENT  OUT  LOOP  800 

C  FOR  TOP  SEMI-INF, LET  YT=3»  AND  COMMENT  OUT  LOOP  860 

C  IF  NOT  USING  SEMI - INF , XL  =  1 »XR  =  X-1 »YT=2 «Y B=Y 

XL  =  1 
XR=X-1 
YB= Y-3 
YT  =  2 

DO  100  I=YT,YB 
DO  2 00  J=XL,XR 
C 

C  EXAMINE  TEMPS  HORIZONTALLY 
R J=R A Y ( I, J, 1 ) 

RJ1=RAY(I,(J+1),1) 

IF((RJ«GT.TISO(l)).AND.(RJl.GT.TISO(l)))  GO  TO  500 
IF((RJ.LT.T1SO(NISO)).AND.(RJ1.LT.TISO(NISO)>>  GO  TO  500 
DO  500  B=l,NISO 

IF< UISO(B) .GT.RJ ) . AND . ( TISO tB > . L T . R J1 ) )  GO  TO  400 
IF<<T1S0CB).LT.RJ).AND.<TIS0(B).GT.RJ1))  GO  TO  400 
IF(TISOIB).EO.RJ)  GO  TO  402 
GO  TO  500 

400  CONTINUE 

401  FORMAT ( IX, »R AY< ♦, 13,13, »,1)=*,F6. 2) 
EXCD=(RAY(I,J,l)-TISO<B))/(RAYCI,J,l)-RAY(I,(J*l),l)) *J 
EXCO=OS*<EXCO-1.) 

COUNT  <B)=COUNT<B)*l 
RECOUNT  IB) 

EXR Y { B, K  >  =EXCO 
EYRY<B»K)=OS*<I-l .) 

CNTPT  =CNTPT  *1 
GO  TO  500 

402  CONTINUE 
COUNT<B)=COUNT<B)*l 
RECOUNT (8) 

EXRY<B,K)=OS*{ J-l  .) 

EYRYCB,K)=0S«<I-1.) 

CNTPT =CNTPT*1 

500  CONTINUE 

C  EXAMINE  TEMPS  VERTICALLY 
RI=RAY«I,J,1> 

RI1=RAYMI-1),J»1> 

IF((RI.GT.TISO(l)).AND.(RIl.GT.TISO(l)))  GO  TO  525 
IFt(RI.LT.TISO(NISO)).AND.«RIl.LT.TISO(NISO)>)  GO  TO  525 
DO  525  B=l,NISO 

1F((TISO(B).GT.RI).ANO.(TISO(B).LT.RI1))  GO  TO  550 
IFUTISOtB).LT.RI).ANO.(TISO(B).GT.RIl)>  GO  TO  550 
GO  TO  525 
550  CONTINUE 

404  F0RMATUX,»RJ=»»F6.2,»  RJ1=«,F6.2) 

EYCD=<RAYUI-1),J,1)-TIS0(B))/CRAY<C  I-1),J,1>- 
6  RAY(I,J,1))*(I-1.) 

EYCO=OS*<EYCO-l.) 

COUNT (B)=COUNT(B)*l 
KsCOUNT(B) 

EYR Y C  B*K  ) sEYCO 
EXRYtB»K)=DS*< J-l.) 

CNTPT=CNTPT*1 
525  CONTINUE 
200  CONTINUE 
100  CONTINUE 

C  800  LOOP  IS  FOR  RIGHT  HAND  SIDE 
00  800  I=VT»YB 
DO  851  Bsl*NlSO 
RI=RAYU,X,1) 

RI1=RAY((I-1),X,1) 

IF((RI.6T.(T|S0(B)>).AND.(RI1 .LT .(TISO(B))))  GO  TO  852 
IF((RI.LT.ITIS0(B))).AND.(RI1.GT.(TIS0(B))>)  GO  TO  852 
1F<RI.EQ.TIS0<8)>  GO  TO  853 


1761) 

1762) 

1763) 

1764) 

1765) 

1766) 

1767) 
1766  ) 

1769) 

1770) 

1771) 

1772) 

1773) 

1774) 

1775) 

1776) 

1777) 

1778) 

1779) 

1780) 

1781) 

1782) 

1783) 

1784) 

1785) 

1786) 

1787) 

1788) 

1789) 

1790) 

1791) 

1792) 

1793) 

1794) 

1795) 

1796) 

1797) 

1798) 

1799) 

1800) 
1801) 
1802) 

1803) 

1804) 

1805) 

1806) 

1807) 

1808) 

1809) 

1810) 
1811) 
1812) 

1813) 

1814) 

1815) 

1816) 

1817) 

1818) 

1819) 

1820) 
1821) 
1822) 

1823) 

1824) 

1825) 

1826) 

1827) 

1828) 

1829) 

1830) 

1831) 

1832) 

1833) 

1834) 

1835) 

1836) 

1837) 

1838) 

1839) 

1840) 

r841> 
842) 
1843) 
1844) 
1845) 
1846) 
1847) 
1848) 
1849) 


852 


853 


851 

800 


96 

94 

92 


91 

99 


95 

90 

666 

93 


LOOP  IS  FOR  TOP  ROW 


862 


863 


861 

860 


60  TO  851 
CONTINUE 

ETCO=U-l.)*tRll-TISO«B)  )/(R  U-R  AY  (  I»X*1 )  ) 
EYC0=DS*(EYCD-1.) 

COUNT (B) =COUNT(B) *1 
K=COUNT(B) 

EYRV(B»K)=EYCO 
EXRY(B*K)=0S*(X-1.) 

CNTPT=CNTPT*1 
60  TO  851 
CONTINUE 

COUNT (B)=C0UNT(B>*1 
K=COUNT(B> 

EYRYCB»K)=DS*(I-1.) 

EXRY(B*K)=DS*(X-1.) 

CNTPT=CNTPT*1 
CONTINUE 
CONTINUE 
860 

XX=X-1 
1=1 

00  860  J=XL*XK 
R J=R A  Y (  1  *  J  *  1 ) 

RJ1=R AYC1.I J*l) ,1 ) 

00  861  B=1*NIS0 

IF((TIS0(B).6T.RJ).AND.(TIS0(B).LT.RJ1)>  60  TO  862 
IF((TIS0(B).LT.RJ>.AND.(TIS0(B).GT.RJ1>)  60  TO  862 
IF<RJ.EO.TISO(B)>  GO  TO  863 
60  TO  861 
CONTINUE 

EXCO=(RAY(I*J*l)-TISO(B))/(RAY(I*J*l)~RAY(I*(J+l)tl))+J 

EXC0=DS»(EXCD-1.) 

COUNT (B)sCOUNT  (B)  *1 
K=COUNT(B) 

EXR  Y (B»K)=EXCO 
EYRY(B*K)=0 
CNTPT=CNTPT*1 
60  TO  861 
CONTINUE 

C0UNT(B)=C0UNT(8)*1 
K=COUNT (B) 

EXRY(B»K)=0S*(J-1.) 

EYRY(B»K)=0 

CNTPT=CNTPT*1 

CONTINUE 

CONTINUE 

T IH=Q  *OELT 

WRITE (13*96)  CNTPT.TIM 

F0RMAT(1X**THE  FOLLOW I NG * *1 5 * ■  POINTS  REPRESENT  TIME=«» 
t  F8 . 4  » •  HOURS:*) 

F0RMAT(1X»F7.2) 

FORMAT (1X*I5) 

DO  666  B=1*NIS0 
L=COUNT (B) 

IF(L.EO.O)  GO  TO  99 
WRITE! 13*91 )  L*B*TISO(B> 

FORMAT (1X*2I6«F6. 2) 

CONTINUE 

00  95  K=1*L 

IF(L.EO.O.)  GO  TO  95 

WRITE (13*90)  (EXR Y(B*K ) *EYRY (B»K ) ) 

CONTINUE 

F0RMAT(1X*2F15.6) 

CONTINUE 

FORMATdX*  *T0TAL  NO.  POINTS  F0UND=**I5) 

RETURN 

END 


AODATA*  DATA-6ATHERIN6  SUBROUTINE 
SUBROUTINE  AODATA 

IMPLICIT  INTEGER*2(A«B»I-L*N»0*W«X*Y*Z) 

IMPLICIT  00U8LE  PRECIS I0N(C-H*M*O*P*R-V) 

C0MN0N/M12I /  X * Y* NISO* Q 

C0NH0N/M12R/  RAY(80(80«3)*TISO(9)*OS*OELT 
COMMON/HI I /  IMAX*A*ITRT«KRNR(4)*ITPC 
C0MH0N/M1R/  01*TDEL*H( 2) * 

A  FLXT(80)*FLXB(80) *FLXL ( 80 ) »FLXR(80)«TMPT(80)» 

8  TMPB(80 ) *TMPL( 80 ) * TMPR( 80 ) *R0PC(2 ) *CPPC (2) »HL (2) *TPC(2) 
AOPOAT  IS  THE  OATA  FILE  FOR  ADIPC 

ALL  UNITS  ARE  IN  METERS«HOURS*CELSIUS 
OELTsTIME  INCREMENT  (HRS) 

IMAXsMAX  NO.  OF  TIME  STEPS  REQUESTED 
QS=DI STANCE  BETWEEN  NOOES  (METERS) 

DIsTHE  MULTIPLE  OF  OS  WHICH  IS  HALF  THE  DISTANCE  TO  INFINIY. 

(DS*DI=HALF  THE  DISTANCE  TO  INFINITY)  DI>1 
TSRFsSURF ACE  TEMP  (OUTSIOE  GRID)  (INFINITE  01ST  AWAY) 
TBTM=BOTTOH  TEMP  (UNDER  6RI0)  (INFINITE  OIST  AWAY) 

TRiTsTEMP  TO  RIGHT  OF  GRID  (INFINITE  OIST) 

TLFtsTEMP  TO  LEFT  OF  GRID  (INFINITE  OIST  AWAY) 

X=  NO.  OF  GRID  NODES  HORIZONTALLY  (X>=5) 

a=  no!  of  §ifM^ntsmXt£rialsl^n  thI*IIid 
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1856) 

1857) 

1858) 

1859) 

1860) 
1861 ) 
1862) 
1865  > 

1864) 

1865) 

1866) 

1867) 

1868) 

1869) 

1870) 

1871) 

1872) 

1873) 

1874) 

1875) 

1876) 

1877) 

1878) 

1879) 

lseo) 

1881) 

1882) 

1863) 

1884) 

1885) 

1886) 

1887) 

1888) 

1889) 

1890) 

1891) 

1892) 

1893) 

1894) 


1895)  C 

1896)  C 

1897)  C 

1898)  C 

1899)  C 


1900) 

1901) 

1902) 

1903) 

1904) 

1905) 

1906) 

1907) 

1908) 

1909) 

1910) 

1911) 

1912) 

1913) 

1914) 

1915) 

1916) 

1917) 

1918) 

1919) 

1920) 

1921) 

1922) 

1923) 

1924) 

1925) 

1926) 

1927) 

1928) 

1929) 

1930) 

1931) 

1932) 

1933) 

1934) 

1935) 

1936) 

1937) 

1938) 

1939) 

1940) 

1941) 

1942) 

1943) 

1944) 

1945) 

1946) 


CP(I,J)=SPECIFIC  HEAT  (W*HR/KG*K> 

TK(  I,J)=THERMAL  CONDUCTIVITY  (W/M*K) 

R0( I»J)=DFNSITY  (KG/M3) 

H(L)=CONVECTION  COEFFICIENT 

CPO  =  $P  EC  I F I C  HEAT  OF  PREVIOUS  TIME  STEP  (FIGURED  IN  MAIN  PROGRAM) 
ROO=DE  NS  I T  Y  OF  PREVIOUS  TIME  STEP  (FIGURED  IN  MAIN  PROGRAM) 
RAY(1»J,1)=PRESENT  NODAL  TEMP 
RAY(1,J»2)=LOCATION  TYPE: 

THIS  PROGRAM  HAS  S EMI -INF  IN  I TE  OPTION  ON  LEFT  « 

BOTTOM,  AND  RIGHT  SIDES  ONLY,  AT  PRESENT. 

ASTERISKS  DENOTE  SITUATIONS  NOT  YET  PROGRAMMED. 

=1  TOP  LEFT  CORNER 
=2  BOTTOM  LEFT  CORNER 
=3  BOTTOM  RIGHT  CORNER 
=4  TOP  RIGHT  CORNER 
=5  INTERIOR  NODE,  FLUCTUATING  TEMP 
=  6  INTERIOR  NODE  ,  CONST  ANT  TEMP 
=7  NODE  ON  CONST  TEMP  BOUNDARY 
=8  CONST  FLUX  BOUNDARY!  TOP 
=27  CONST  FLUX  BOUNDARY:  LEFT 
=26  CONST  FLUX  BOUNDARY:  BOTTOM 
=29  CONST  FLUX  BOUNDARY:  RIGHT 
=9  CONVECTIVE  BOUNDARY:  TOP 
=30  CONVECTIVE  BOUNDARY:  LEFT 
=31  CONVECTIVE  BOUNDARY:  BOTTOM 
=32  CONVECTIVE  BOUNDARY!  RIGHT 
=1C  NODE  ON  SEMI -INFINITE  BOUNDARY 

NODE  ADUACENT  TO  LEFT  SEMI-INF  BNDRY 
NODE  ADUACENT  TO  BOTTOM  SEMI-INF  BNDARY 
NODE  ADJACENT  TO  RIGHT  SEMI-INF  BNDRY 
NODE  ADJACENT  TO  TOP  SEMI-INF  BNDRY 


=  11 
=  12 
=  13 
=  14 


THE  NEXT  FIVE  CONDITIONS  NEEDED  WITH  SEMI-INF  CORNER 


:15 
=  16 

=  17 
=  18 
=  19 


=  20 
=  21 
=  22 
=23 


SQUARE  NODE  ADJACENT  TO  TWO  SEMI-INF  SIDES 
SEMI-INF  NODE  ABOVE  SEM.-INF  CORNER  NODE 
(FOR  BOTTOM  LEFT  &  RIGHT  ONLY,  AT  PRESENT) 

SEMI-INF  NODE  TO  LEFT  OF  BOTTOM  RIGHT  SEMI-INF  CORNER  NODE 
SEMI-INF  NODE  BELOW  SEMI-INF  CORNER  NODE  * 

SEMI-INF  NODE  TO  RIGHT  OF  BOTTOM  LEFT  SEMI-INF  CORNER  NODE 
THE  NEXT  FOUR  CONDTNS  NEEDED  WITH  CORNER  W  ONE  SIDE  SEMI-INF, 

ONE  SIDE  CONST  FLUX 

TOP  CONST  FLUX  NODE  ADJACENT  TO  RIGHT  SEMI-INF 
RIGHT  SIDE  CONST  FLUX  NODE  ADJACENT  TO  BOTTOM  SEMI-INF 
TOP  CONST  FLUX  NODE  ADJACENT  TO  LEFT  SEMI-INF 
LEFT  SIDE  CONST  FLUX  NODE  ADJACENT  TO  BOTTOM  SEMI-INF 
THE  NEXT  THREE  CONDTNS  NEEDED  WITH  CORNER  W  ONE  SIDE  SEMI-INF, 

ONE  SIDE  CONVECTIVE 

=24  RIGHT  SIDE  CONVECT  NODE  ADJ  TO  BOTTOM  SEMI-INF. 

=25  TOP  CONVECT  NOOE  ADJ  TO  LEFT  SEMI-INF* 

LEET  SIDE  CONST  FLUX  NODE  ADJACENT  TO  BOTTOM  SEMI-INF* 
RAY(I«J«3)=  INDEX  OF  MATERIAL  (RANGES  FROM  1  TO  A) 

KRNR(1)=T OP  LEFT  CORNER 

KRNR(2)=B0TT0M  LEFT  CORNER 

KRNR ( 3 ) =BOTTOM  RIGHT  CORNER 

KRNR(4)=T0P  RIGHT  CORNER 

KRNR ( L ) =1  CONSTANT  TEMP  CORNER  NODE 

=2  CONSTANT  FLUX  (ON  BOTH  S1DEJ)  CORNER  NODE 
=3  CONVECTIVE  (ON  BOTH  SIDES)  CORNER  NODE 
=4  SEMI-INFINITE  (ON  BOTH  SIDES)  CORNER  NODE 
=  5  VERTICAL  S I DE  =  CONVECT ,  HORIZ  =  CONST  FLUX 
=6  VERTICAL  SIDE=CONST  FLUX,  HORIZ  =CONVECT 
=7  VERTICAL  SIOE  =CONST  FLUX,  HORIZ  =SEMI-INF 
=8  VERTICAL  S I 0E=SEM1 -I NF ,  HORIZ=CONST  FLUX 
(TOP  LEFT  &  TOP  RIGHT  ONLY  AT  PRESENT) 

=9  VERTICAL  SIDE  =CONVECT*  HOR I Z=SEMI -I NF  * 

=1°  VERTICAL  SIOE=SEMI-INF,  HORIZ  =  CONVECT* 

TEMPI  1 , J) =  TEMP  AT  TIME  ( T *DELT ) 

ISTAT  (I  *J)  =  INDEX  OF  STATE  FOR  TEMPO, J) 

=1  FROZEN 
=2  UNFROZEN 

=3  UNDERGOING  PHASE  CHANGE 
FLX  T ( J) =FLU  X  FROM  TOF  OF  GRID 
FLXB(J)=FLUX  FROM  BOTTOM  OF  GRID 
FLXL(I)=FLUX  FROM  LEFT  SIDE  OF  GRID 
FLXR(I)=FLUX  FROM  RIGHT  SIDE  OF  GRID 
TMP T ( J ) =TEMP  OUTSIDE  GRID  FOR  TOP 
THPB ( J ) =TEMP  OUTSIDE  GRID  FOR  BOTTOM 
TMPL(I)=TEMP  OUTSIDE  GRID  FOR  LEFT 
TNPR(I)=TEMP  OUTSIDE  GRID  FOR  RIGHT 
N ISO=  NO.  OF  ISOTHERMS  TO  BE  PLOTTED 
TMAX=  TEMP  OF  HOTTEST  ISOTHERM 
TMIN=  TEMP  OF  COLDEST  ISOTHERM 
TISO(B>=  ARRAY  CONTAINING  ISOTHERM  TEMPS 
TISO(l)  HAS  THE  HOTTEST  ISOTHERM 
COUNT(B)=COUNTER  FOR  NO.ELTS  IN  EACH  ISOTHERM 
ITRT=NO.  TIME  STEPS  BEFORE  TEMPS  PRINTED 
I  TPC  =  NO .  TIME  STEPS  BEFORE  ISOTHM  IS  CALLED 
FOR  THE  NEXT  4  STM TS »K =R AY < I , J , 3) 

ROPC(K)=OENSITY  AT  PHASE  CHANGE  FRONT (KG/M3) 

HL(K)=LATENT  HEAT  OF  FUSION  (W*HR/KG*K) 

TPC(K)=TEMPERATURE  OF  PHASE  CHANGE  (»C> 

CPPC(K)=SPECIFIC  HEAT  DURING  PHASE  CHANGE 
TOEL  =  TOT AL  TEMP  RANGE  FOR  TPC 

(EX!  IF  TPC=0»»TOEL=2,PHSCH  OCCURS  FROM  -1»  TO  !•) 
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1993) 

1994) 

1995) 

1996) 

1997) 

1998) 

1999) 

2000) 
2001) 
2002) 

2003) 

2004) 

2005) 

2006) 

2007) 

2008) 

2009) 

2010) 

18l2> 
2013) 
2  014) 

2015) 

2016) 

2017) 

2018) 

2019) 

2020) 
2021) 
2022) 

2023) 

2024) 

2025) 

nm 

2028) 

ISIS! 


C  INITIALIZE  VARIABLES  ANO  ARRAYS 
TSRF=-4. 67000 
TBTM=4.6700 
TLFT=4. 67000 
TRIT=4. 67000 
DS=.0050D0 
DELT=. 00250D0 
01=50.000 
A=1 

INAX=1200 
X=3 
Y  =25 

TDEL=1.000 
NIS0=1 
I TRT=40 
I TPC=10 
C 

C  00  NOT  CHANGE  THE  FOLLOWING  6  LINES. 

ETIME=ITRT*DELT 

MTIME=IMAX*DELT 

WRITE (1.3)  ITRT  .ETIME.IMAX.MT1ME 

3  F OR HA T <lXt ‘TEMPER ATURES  WILL  BE  PRINTEO  EVERY*. 13. •  ITNS*./. 

C  IX.*  (EVERY*. F12.4**  HRS.)*«/*1X* 


C  *M AX  NO.  I TNS=  * » 1 4 . • 


<=*«F12.4»*  HRS)*) 


SET  UP  MATERIAL  PROPERTIES  FOR  EACH  MATERIAL. 

(YOU  SET  THE  VALUES  OF  K  A  CP  &  RO  IN  MAIN  PROGRAM.) 

H(A)=O.DO 

PROPERTIES  FOR  PHASE  CHANGE: 

ROPC(1>=917.000 
HL( 1 >=93.000 
TPC (1 ) =0 • OOO 

THE  FOLLOWING  SPECIFIC  HEAT  SHOULD  INCLUDE  LATENT  HEAT. 

YOU  CALCULATE  IT  BY:  CP= ( CP (FR OZ )*CP( UNFROZ) )/2  ♦  HL(K)/T0EL 
00  2  K=1 « A 
CPPC(K)=093. 87000 
CONTINUE 

SPECIFY  THE  ISOTHERMS  TO  BE  LOCATED  IN  SUBROUTINE  ISOTHM 
IN  ORDER.  WITH  TISO(l)  THE  HOTTEST. 

T ISO( 1 )=0. 0D0 

INDICATE  THE  FLUXES  FROM  THE  SIDES  OF  THE  GRID. 

THIS  IS  USED  ONLY  FOR  THE  CONSTANT  FLUX  BOUND .CONDITION. 

DO  11  J=1.X 
FLXT ( J)=0«00 
FLXB( J>=0.D0 
11  CONTINUE 
00  12  1=1. Y 
FLXL( I )=0.00 
FLXR ( I )=0 .00 

12  CONTINUE 
C 

C  INDICATE  THE  TEMPERATURES  OUTSIDE  THE  GRID. 

C  USED  FOR  SEMI-INF  &  CONVECT  BOUNDARYS 

DO  13  J=1«X 
TMPT ( J ) =TSRF 
TMPB( J)=TBTM 

13  CONTINUE 

DO  14  1=1. Y 
TMPL(I)=TLFT 
TMPR( I)=TRIT 

14  CONTINUE 


SET  UP  RAY(I.J»3>  INOEX  OF  MATERIALS 

THERE  ARE  A  MATERIALS.  LET  THE  SURROUNDING 
MATERIAL  BE  THE  *ATH*  MATERIAL.  START 
WITH  *1*. 

THIS  LOOP  ASSIGNS  THE  • ATH •  MATERIAL  TO 
THE  REMAINING  NODES. 

00  10  J=1 • X 


20 

10 


J=l») 
DO  20  1=1. Y 
RAY(1«J.3)=A 
CONTINUE 
CONTINUE 


C  SET  UP  RAV( 1 . J.2)  NODAL  LOCATION  TYPE 
C  TOP  BOUNDARY 
DO  40  J=1 .X 
R AY(1 »J.2) =7 
40  CONTINUE 
C  BOTTOM  BOUNDARY 
DO  50  J=l.X 
RAY(Y*J.2)=10 
SO  CONTINUE 
C  LEFT  BOUNDARY 
DO  60  1=1. Y 
RAY(I*1.2>=27 
60  CONTINUE 
C  RIOHT  BOUNOARY 
00  70  I =1. Y 
R AY ( 1 *X .2) =29 
70  CONTINUE 

C  SET  UP  CORNER  CONDITIONS 
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✓  / 


2  041) 

2042) 

2043) 

2044) 

2045) 

2046) 

2047) 

2048) 

2049) 

2050) 

2051) 

2052) 

2053) 

2054) 

2055) 

2056) 

2057) 

2058) 

2059) 

2060) 
2061  ) 
2062) 

2063) 

2064) 

2065) 

2066) 

2067) 

2068) 

2069) 

2070) 

2071) 

2072) 

2073) 

2074) 

2075) 

2076) 

2077) 

2078) 

2079) 

2080) 
2081) 
2082) 

2083) 

2084) 

2085) 

2086) 

2087) 

2088) 

2089) 

2090) 

2091) 

2092) 

2093) 

2094) 

2095) 

2096) 

2097) 
2098> 

2099) 

2100) 
2101) 
2102) 

2103) 

2104) 

2105) 

2106) 

2107) 

2108) 

2109) 

2110) 
2111) 
2112) 

2113) 

2114) 

2115) 

2116) 

2117) 

2118) 

2119) 

2120) 
2121) 
2122) 

2123) 

2124) 

2125) 

2126) 

2127) 

2128) 

2129) 

2130) 

2131) 

2132) 

2133) 

2134) 


--DON’T  CHANGE  THE  "RA Y( I .U.2) *  STATEMENTS, 
RAY(1.1«2>=1 
R AY( Y« 1 .2) =2 
RAY(Y,X,2>=3 
RAY(1.X.2>=4 

TOP  LEFT  CORNER 
KRNR C 1 ) = 1 
RAY(1,1.1)=TSRF 

BOTTOM  LEFT  CORNER 
KRNR ( 2  >  =7 
RAY(Y.l.l) =TBTM 

BOTTOM  RIGHT  CORNER 
KRNR  <  3>=7 
RAY(Y,X,1)=TBTM 

TOP  RIGHT  CORNER 
KRNR  t  4>  =  1 
RAY(1,X.1)=TSRF 

:  INTERIOR  NODES  GIVEN  FLUCTUATING  TEMP 
:  DO  NOT  CHANGE  THE  80.90  LOOP. 

YY=Y-1 

XX=X-1 

DO  8C  J=2.XX 
DO  90  I  =2 » Y Y 
RAY(I«J«2)=5 
90  CONTINUE 
80  CONTINUE 


CONSTANT  INTERIOR  NODES. 

INDICATE  R AY ( I » J. 2 ) =6  IF  RELEVANT 

IF  RELEVENT.  SPECIFY  SPECIAL  RAYCI,J,2)  CONDITIONS 
TO  BE  USED  WITH  SEMI-INF  BOUNDARY. 

XX=X-1 

YY=Y-1 

DO  210  J=2.XX 
R AY(YY«J.2>=12 
210  CONTINUE 
:  RAV(YY.2«2)=15 

RAY(YY.1.2)=23 
R AY ( Y.2.2) =19 
RAY(YY.XX.2)=15 
RAY ( VY.X.2) =21 


SET  UP  RAV(I.J.l)  PRESENT  NODAL  TEMP 
DO  NOT  ALTER  THIS 
XX=X-1 
YY=Y-1 

DO  110  J=2.XX 
RAY(1»J»1)=TSRF 
RAY(Y,J,1)=TBTM 
110  CONTINUE 

DO  120  1=2. YY 
RAYU.1.1)=TLFT 
RAY(I,X.1)=TRIT 
120  CONTINUE 
XX=X-1 
Y Y=Y-1 

00  130  J  =  2  *  XX 
DO  140  1=2. YY 

C  COMMENT  OUT  THE  INITIAL  TEMP  DIST  THAT  YOU  DON’T  WANT: 
C  NOOES  ASSIGNED  LINEAR  TEMP  DIST  VERTICALLY 
C  RAY(I.J»1)=TSRF-((TSRF-T8TM)/(Y-1))»(I-1) 

C  NODES  ASSIGNED  UNIFORM  TEMP  (INDICATE  THE  TEMP:) 
RAY(I.J«1)=4.670D0 
140  CONTINUE 
130  CONTINUE 


INTERIOR  NODES 

INDICATE  THE  TEMPERATURES  FOR  RAY(l.J.l) 

IF  THEY  ARE  DIFFERENT  THAN  ASSIGNED  IN  THE  130.140 


DO  NOT  CHANGE  THE  FOLLOWING  17  STATEMENTS. 

WR I TE (5. 131 )  DS.DELT.OI. TDEL 

131  FORMATC1X.4F10.5) 

WRITE (5 .132)  A.X» Y, NISO. ITRT.I MAX. 1TPC 

132  FORMAT ( IX. 715) 

WRITE (5* 136)  ((ROPC(K) .CPPC ( K ) .HL ( K ) . TPC (K ) ) .K =1 . A) 
136  F0RMAT(1X.4F10.5) 

WRITE(5.134)(T1S0(B).6=1.NIS0) 

134  F0RHAT(1X,F7.2) 

WRITE (5.133)  (H(L)«L=1.A) 

133  FORMAT(1X*F10.5) 

WRITE (5. 151)  ((RAY(I.J«3).J=1«X).I=1»Y) 

WRITE (5. 151)  ((RAY(I»U«2).J=1»X),I=1.Y> 

WRITE (5* 151)  ((RAY(I.J.1).J=1.X)*1=1*Y) 

WRITE (5. 152 )  (KRNR(J)«J=1«4> 

152  F0RHAT(1X.4I2> 

151  FORMAT! 1X.17F7. 2) 

WRITE (5.151)  ( FLX T  (J).U  =  1*X) 

WRITE (5. 151)  (FLXB(U). J=1,X) 

WRITE (5. 151)  (FLXL ( I).I=1.V) 

WRITE (5. 151)  (FLXR(I).I=1.Y) 
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FUe  ADPOUT 

(Lists  the  initial  data  in  readable  form  and  the 
final  values  from  the  run  [printed  by  ADIPCJ.) 

DATA  FOR  THIS  RUN  OF  AOIPC: 


OS 

0.00500 


0.00250  50.00000 


T  DEL 

1.00000 


Y  NISO  ITRT  I H  A  X  ITPC 
25  I  40  1200  10 


ROPC  <  K  >  CPPC  < K  )  HL  <  K  >  TPCOO 

917.00000  93.87000  93.00000  0.00000 

TIS0(fi) ,B=ltNIS0: 

0.00 

H  ( K  > 

0.00030 

KRNR<1>=  1  KRNR<2>=  7  KRNR<3)=  7  KRNR(4>=  1 

FL  XT ( J )  FLXfltJ)  TMPT(J)  TMPBtJ) 

0.00  0.00  -4.67  4.67 

0.00  0.00  -4.b7  4.67 

0.00  0.00  -4.67  4.67 

FL  XL  1 1  >  FLXR(I)  TMPLtI)  TMPR(I> 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

O.uO  0.00  4.67  4.67 

C.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 

0.00  0.00  4.67  4.67 


3 1  J ) 

TMPTtJ) 

TMPBtJ) 

0.0  0 

-4.67 

4.67 

0.00 

-4.67 

4.67 

0.00 

-4.67 

4.67 

til) 

T  MPL (  I  ) 

TMPR  (  I ) 

0.0  0 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

0.00 

4.67 

4.67 

RAY(I»Jtl) 

>  TEMPERAT 

-4.67  -4.67 

-4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.6T 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

4.67 

R AYt It J*2) 

NOOAL  LOI 

1.00 

7.00 

4.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

62 


27.00 

5.0  0 

29.00 

27.00 

5 . 0  C 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29.00 

27.00 

5.00 

29. CO 

23.00 

12.00 

21.00 

2.00 

10.03 

3.00 

RAYU.J.3)  NODAL  MATERIAL  TYPE 


1.00 

1.00 

1.00 

1.00 

1.03 

1.00 

1.00 

1 .  G  0 

1.00 

1.00 

1.03 

1.00 

1  .00 

1.03 

1.00 

1.00 

1.00 

1.03 

1.00 

1.30 

1 .  P  0 

1.00 

1.00 

1 . 0  o 

1.00 

1.03 

1.00 

l.CO 

1.03 

1.00 

1.00 

1.0G 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.03 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1  •  u  0 

1.00 

l.CO 

1.00 

1.00 

1.00 

1.00 

1.03 

1.00 

TEMPERATURES  WILL  BE  PRINTED  EVERY  AO  TIME  STEPS. 
ISOTHERMS  WILL  BE  LOCATED  EVERY  10  TIME  STEPS. 

BOUNDARY  CONDITIONS: 

TOP  BOUNDARY  CONSTANT  TEMPERATURE 

TOP  LEFT  CORNER  CONSTANT  TEMPERATURE 

LEFT  BOUNDARY  CONSTANT  HEAT  FLUX 

BOTTOM  LEFT  CORNER  VERT  CONST  FLUX — HORIZ  SEMI-INF 
BOTTOM  BOUNDARY  SEMI -INF INITE 

BOTTOM  RIGHT  CORNERVERT  CONST  FLUX--HORIZ  SEMI-INF 

RIGHT  BOUNDARY  CONSTANT  HEAT  FLUX 

TOP  RIGHT  CORNER  CONSTANT  TEMPERATURE 

NUMBER  OF  TIME  STEPS=  1200 


TEMPERATURES  AFTER1200  TIME  STEPS  <  3.000  HRS): 


A  .67 

-A. 67 

-4.67 

3. 7A 

-3.7  A 

-3.74 

2.82 

-2.82 

-2.82 

1.89 

-1.89 

-1.89 

0.97 

-0.97 

-0.97 

0  •  0  A 

-0.0  A 

-  0 . 0  A 

0.38 

0.38 

0.38 

0.65 

0.65 

0.65 

1.07 

1.07 

1.07 

1.A6 

1.46 

1.46 

1 .84 

1  .84 

1.84 

2.18 

2.18 

2.18 

2.50 

2.50 

2.50 

2.79 

2.79 

2.79 

3 . 0  A 

3 . 0  A 

3. 04 

3.27 

3.27 

3.27 

3.A7 

3.47 

3.47 

3 . 6  A 

3.64 

3.64 

3.79 

3.79 

3.79 

3.91 

3.91 

3.91 

4.00 

4.00 

4.00 

4.07 

4.07 

4.07 

4.13 

4.13 

4.  13 

4.16 

4.16 

4.16 

*  A 


FINAL  TK<I»J): 


2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

2.21 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.56 

0.58 

0.58 

0.58 

0.58 

0.58 

0.56 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.56 

0.58 

0.58 

0.56 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.56 

*AL  CPCI.J): 

0.58 

0.58 

0.56 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.58 

0.56 

0.58 

0.58 

0.58 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

1.16 

FINAL  RO(I.J): 

917.00  917.00  917.00 
917.00  947.00  917.00 
917.00  917.00  917.00 
IH*82  917.00  917.00 
917.00  917.00  917.00 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 
998«20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998120 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

998.20  998.20  998.20 

final  istatu.j): 

1  1  1 
1  1  1 
1  1  1 
1  1  1 
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Sample  of  isotherm  locations  output  from  ADIPC  for  the 
Neuman  problem  (file  POINT  1) 


THE  FOLLOWING  3 

3  1  0.00 

0. 000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 

THE  FOLLOWING  3 

3  1  0.00 

0.000000 
0.005000 
0.010000 


POINTS  REPRESENT  TINE: 

0.00463a 

0.004632 

0.00463a 

POINTS  REPRESENT  TINE: 

0*004727 

0.004727 

0.004727 

POINTS  REPRESENT  TINE: 

0.004830 

0.004830 

0.004830 

POINTS  REPRESENT  TINE: 

0.004940 

0.004940 

0.004940 

POINTS  REPRESENT  TINE: 

0.005197 

0.005197 

0.005197 

POINTS  REPRESENT  TINE: 

0.005661 

0.005661 

0.005661 

POINTS  REPRESENT  TINE: 

0.006161 

0.006161 

0.006161 

POINTS  REPRESENT  TINE= 

0.006685 

0.006685 

0.006685 

POINTS  REPRESENT  TINE: 

0.007223 

0.007223 

0.007223 

POINTS  REPRESENT  TINE: 

0.009329 

0.009329 

0.009329 

POINTS  REPRESENT  TINE: 

0.009413 

0.009413 

0.009413 

POINTS  REPRESENT  TINE: 

0.009514 

0.009514 

0.009514 

POINTS  REPRESENT  TINE: 

0.009618 

0.009618 

0.009618 

POINTS  REPRESENT  TINE: 

0.009724 

0.009724 

0.009724 

POINTS  REPRESENT  TINE: 

0.009831 

0.009831 

0.009831 

POINTS  REPRESENT  TINE: 

0.009941 

0.009941 

0.009941 

POINTS  REPRESENT  TINE: 

0.010158 

0.010158 

0.010158 


0.0250  HOURS: 

0.0500  hours: 

0.0750  hours: 

o.iooo  hours: 

0.1250  hours: 

0.1500  hours: 

0.1750  hours: 

o.2ooo  hours: 

0.2250  hours: 

0.2500  hours: 

0.2750  hours: 

0.3000  hours: 

0.3250  hours: 

0.3500  hours: 

0.3750  hours: 

0.4000  hours: 

0.4250  hours: 


A  facsimile  catalog  card  in  Library  of  Congress  MARC 
format  is  reproduced  below. 
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